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We present a new numerical relativity code designed for simulations of compact binaries involving 
matter. The code is an upgrade of the BAM code to include general relativistic hydrodynamics and 
implements state-of-the-art high-resolution-shock-capturing schemes on a hierarchy of mesh refined 
Cartesian grids with moving boxes. We test and validate the code in a series of standard experiments 
involving single neutron star spacetimes. We present test evolutions of quasi-equilibrium equal-mass 
irrotational binary neutron star configurations in quasi-circular orbits which describe the late inspiral 
to merger phases. Neutron star matter is modeled as a zero-temperature fluid; thermal effects can 
be included by means of a simple ideal-gas prescription. We analyze the impact that the use of 
different values of damping parameter in the Gamma-driver shift condition has on the dynamics 
of the system. The use of different reconstruction schemes and their impact in the post-merger 
dynamics is investigated. We compute and characterize the gravitational radiation emitted by the 
system. Self-convergence of the waves is tested, and we consistently estimate error-bars on the 
numerically generated waveforms in the inspiral phase. 
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I. INTRODUCTION 

Binary neutron stars (BNS) are among the most 
promising sources of gravitational waves (GWs) for 
ground-based interferometers of present and future gen- 
eration and are also at the origin of the powerful 
electromagnetic astrophysical phenomena, in particular 
short-gamma-ray- bursts (SGRBs) HQ. While other as- 
pects of BNS (and neutron star, NS) physics are cer- 
tainly interesting, these two topics represent to date one 
of the most exciting observational and theoretical chal- 
lenge. On the one hand, the detection of GWs emitted 
during the last orbits of a merger process is expected to 
convey unique information about the nature of matter at 
supra-nuclear densities which is largely unknown, e.g. Q. 
On the other hand SGRBs are ultra-relativistic outflows 
likely to be injected during the post-merger dynamics 0- 
but neither a self-consistent model nor a simulation 
are yet available to provide the conditions for the "central 
engine" H3, E3 - 



A complete theoretical modeling of the late inspiral 
and merger phase is possible only by means of numerical 
relativity (NR) simulations. While BNS simulations have 
a longer history, see e.g. [l2| - [l3 | for recent reviews, the 
first NR simulation was performed in fl5l ]. some years 
before the first complete simulations of coalescing binary 
black- holes (BBHs) [l6Hl8|. Nowadays a number of NR 
groups are performing BNS simulations [l9l - [25j . 

Recent work investigated the evolution of irrotational, 
circularized, e qual and unequal mass binaries without 
magnetic fields |19l . l20l I26rl28j . In these works the dynam- 
ics of the system, mainly dependent on the initial masses 
of the stars involved, is explored in detail with particu- 
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lar attention on the final product of the merger: either 
a hyper-massivc-NS (HMNS) (which eventually collapses 
on dynamical timescales) or a black-hole (BH) with an 
accretion disk configuration resulting from a prompt col- 
lapse. The gravitational radiation emitted by the systems 
has been characterized. 

Electromagnetic fields in NR simulations of BNS have 
been considered in @, HH, [H [H Ho[ . Their impact on 
GWs during the inspiral has been found negligible for as- 
trophysically motivated intensities [H,[29|,[3(|, while they 
have certainly a major role in the post-merger phase for 
several astrophysical scenarios like SGRBs. Up to now 
electromagnetic fields have been considered in full gen- 
eral relativity (GR) simulations only within the frame- 
work of ideal general relativistic magneto- hydrodynamics 
(GRMHD), i.e. in the limit of infinite electrical conduc- 
tivity. Some efforts towards resistive GRMHD are ongo- 
ing |31| - 

In all cases mentioned above the studies are based 
on the numerical solution of the ideal general relativis- 
tic hydrodynamics equations (GRHD), or GRMHD, cou- 
pled to some 3+1 hyperbolic formulations (BSSNOK (32r - 
□3 or GHG (25|) of GR. Most of the NR results are 
obtained with simplified treatments of the NS interior, 
namely ideal gas, polytropic or piccewise polytropic equa- 
tions of state (EoS). In [36[ (and following works) zero- 
temperature (cold) "realistic" EoS are also employed. 
Thermal effects are expected to be relevant in the post- 
merger phase. They are taken into account in an approx- 
imate way using hybrid EoS [3H [37j , or in simulations 
based on approximations of GR which instead focus more 
on microphysical aspects (see e.g. [38M40j ) . The model- 
ing of microphysics is particularly important in combi- 
nation electromagnetic fields in order to model SGRBs. 
Transport phenomena and neutrino physics [41H43| are 
currently not included in full GR BNS simulations (see 
however the very recent work in [44^. 

We also mention that mixed binaries, i.e. binary sys- 
tem composed of a black-hole and a neutron star, arc 
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currently under investi gat ion with the same techniques 
employed for BNSs |45l - [53 |. 

While some of the above mentioned aspects of BNS 
physics are quite well understood, the complexity of the 
physical and technical problem certainly requires more 
investigations. On the one hand there is the need of im- 
proving the microphysical modeling of the post-merger 
phase, on the other hand there is the necessity to inves- 
tigate the initial configuration space and to improve the 
accuracy of the waveforms, cf. [l3( . 

The latter point is particularly important for future 
GW detection, see for example [53J, also in view of a 
starting a program of collaboration between NR and an- 
alytical relativity community @, [53 - f57j . In contrast to 
BBH simulations, where a precise waveform analysis is 
routinely performed (see e.g. [H, HH), the accuracy of 
NR waveform from BNS has been poorly investigated so 
far. Ref. [(3(| reports a first error analysis of BNS wave- 
forms. The curvature (^4) waveforms are found to con- 
verge at order 1.8 in the inspiral phase and at 1.2 after 
the merger. The waveforms are aligned to perform the 
convergence test [6ll ] . An error budget regarding the grid 
setting is also presented. In [H?} a similar error budget 
regarding finite extraction, thermal and resolution effects 
is presented. An estimate of the phase and amplitude er- 
rors during inspiral is given based on a convergence rate 
assumed from other simulations [6(| . The main source of 
error has been found to be the finite resolution. We are 
not aware of other detailed investigations regarding con- 
vergence and precise error estimates on the numerically 
generated GWs. 

In this paper we present a new code to simulate non- 
vacuum spacetimes in full GR. The code is an extension 
of the BAM code developed at Jena and elsewhere (62l — 
HH for numerical studies of multiple black holes space- 
times with adaptive mesh refinement techniques. The 
BAM code has been upgraded in order to solve the 
flux-conservative Eulerian formulation of ideal GRHD 
equations [65[ coupled to the Einstein system. We de- 
scribe in detail the equations and the implementation of 
HRSC scheme in BAM. The code allows the use of hy- 
brid EoS composed of a cold part, generically provided 
by tables, and a thermal part modeled with an ideal gas 
EoS [1^, [37} . A simple thermodynamical consistent pro- 
cedure for table interpolation is employed. 

We validate the code against a number of stringent 
tests involving single-star spacetimes. Each test permits 
the verification of a part of the numerical algorithm. 
Convergence and constraint violation in particular are 
discussed in some detail considering different reconstruc- 
tion methods. 

We present our first results concerning the simulation 
of gravitational radiation emitted from BNS evolutions. 
Since we do not simulate magnetic fields, our main in- 
terest is related to the GWs emitted during the inspiral 
phase. We focus, as a test case, on a binary already 
considered in the literature [22|, H3, |66| . The initial ir- 
rotational configuration in quasi-circular orbit has been 



evolved for about three orbits, the merger and the post- 
merger phase with both a cold and a hot EoS. We in- 
vestigate the effect of the use of different reconstruc- 
tion methods on the dynamics. We report new results 
concerning the role of the damping parameter, 77, in the 
Gamma-driver shift condition with regard to the dynam- 
ics. Waveforms are extracted from the simulation via a 
standard algorithm based on the Newman-Penrosc scalar 
-04. The actual GWs (metric waveforms) are additionally 
reconstructed with two different methods, namely a time- 
domain [63, [68[ and a spectral [6|| integration. The two 
methods are compared. We discuss the convergence of 
the numerically extracted GWs. Using different conver- 
gence series we provide the first consistent estimate of 
the errors on phase and amplitude in BNS simulations. 
This work is our first contribution to the study of BNS 
spacetime by NR simulations; the method presented here 
will serve as a future reference. 

The structure of the paper is as follows. In Sec. UT1 we 
review the basic equations solved in BAM. In Sec. IIIII 
we review the description of the NS matter within BNS 
simulations in NR. In Sec. IIVI we summarize the numer- 
ical method adopted for solving the GRHD equations as 
well as the singularity and vacuum treatment in BAM. 
In Sec. |V] we present our results concerning single star 
spacetimes. In Sec. I VII we present our results on simula- 
tions of BNSs. 

Dimcnsionless units c = G = M@ = 1 are employed. 
Times and lengths are often expressed in ms and km 
to facilitate comparison with the literature. Indexes 
a, 6, c, ... run from to 3, indexes i,j,k,... from 1 to 3. 

II. EQUATIONS 

In this section we review the equations solved in BAM. 
We assume the usual 3+1 decomposition of spacetime; 
the metric reads 

ds 2 = -{a 2 - ft^dt 2 + 2p i dx i dt + ^jdx'x , (1) 

where a and ft 1 are the lapse and shift vector and 7,-j the 
spatial 3-metric. The Einstein equations are formulated 
in the strongly hyperbolic BSSNOK form and presented 
in Sec. ITTAl While the BSSNOK system is described in 
several textbooks we describe it again in the following 
for completeness. The equation are explicitly written for 
the x _ BSSNOK system not given in [621 ] ; we point out 
some relevant detail for the implementation and correct 
some minor typos. GRHD equations are given in Sec. Ill Bl 
following the flux-conservative formulation of 65] . 

A. Metric 

The BSSNOK formalism assumes the conformal de- 
composition of the 3-metric, 

7*3 = e _4 *7ij and <t> = T77 mdet 7« > ( 2 ) 
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where 7^ is the conformal 3-metric and <f> the conformal 
factor. Following [l8[ we introduce the variable 

X = e- i4> . (3) 

The extrinsic curvature K~ij of the spatial hypcrsurfaccs 
is decomposed as 

An =X[K ij - I^k) and K = 7 «JT ij . (4) 



(dt - Cp) X 
(d t - C f3 )j tj 
(d t - Cp)Aij 

(d t -£p)K 
dt? 1 



Above Cp is the Lie derivative along the shift vector, 
Di the covariant derivative associated with 7^, R^ is 
the Ricci tensor, see Eq. (18)-(20) in (62j , and the York- 
ADM quantities are defined in Sec. Ill Bl The terms pro- 
portional to padm in Eqs. (0 and (JH| appear because 
the trace of Rij is obtained by substituting the Hamilto- 
nian constraint as usual, and in Eq. (J7J we choose to use 
the vacuum constraint in (i? ij ) TF and to add the missing 
term separately. The variables 

r = 7> fc f;. fc , (10) 

defined in term of the Christoffcl symbols of the confor- 
mal metric, are promoted to new evolution variables. 
The constraints are 

d ee %fi - f k d k % , (11) 
H ee R- A^A l3 + \k 2 - 167TPADM , (12) 

AP ee Bjtfl + f* jk A* - p ij djK - logx • 

(13) 

The gauge is specified by the 1+log lapse and Gamma- 
driver-shift [nHzi: 

d t a = P l d a t - a 2 p L K , (14) 
9 t P = Hsf i -TlP + P i d i p . (15) 

The parameters are fixed to ps = 3/4 and 77 = 0.3 if not 
stated otherwise. 



The evolution system reads 



(5) 
(6) 
(7) 

(8) 

(9) 



B. Matter 

We assume matter composed of a single particle species 
(simple fluid) and described by the perfect fluid stress- 
energy tensor 

T ab = phu a u b + pg ab , (16) 

where p is the rest-mass density, e is the specific internal 
energy, h = 1 + e + p/pis the specific enthalpy, p is the 
pressure, and u a is the 4- velocity (u a u a = —1) of the 
fluid. The total energy density is given by e = p(l + e). 

The GRHD equations for the perfect fluid matter 
(ideal GRHD) are the local conservation law for the 
energy-momentum tensor, the conservation law for the 
baryon number, and the equation of state of the fluid: 

V a T ab = , (17) 
V a (pu a ) = , (18) 
P(p,e) = p . (19) 

Following [65| we rewrite Eqs. ([T7|) - ([T5| in first-order 
flux-conservative form, 

d t q + d i f^(q) = s(q) , (20) 

by introducing the conservative variables 

q = Vl{D, S k ,r}, (21) 
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where 

D = Wp , 

S k = W 2 phv k , (22) 

r = (W 2 ph -p) - D . 

The simple physical interpretation of these variables is 
that they represent the rest-mass density (D), the mo- 
mentum density (Sk) and an internal energy (r = padm— 



D) as viewed by Eulerian observers. Above v % is the fluid 
velocity measured by the Eulerian observer with 

W a a \u u J 

W is the Lorentz factor between the fluid frame and the 
Eulerian observer, W = l/y/l — v 2 , with v 2 = jijv*^ . 
The fluxes in Eq. (J20J) are 



fi i > = y/=9{D[v i - 



,S k [v l - 



+ p5\ ,t «*- 



0* 



■ pv 



while the source terms are 

s = V~g{^T ah (d a9ak - r s ab g Sk ) ,a (T a0 d a ]na - T ab T° ab )} 
= J—g |o,T 00 Q^ATtj ~ <*0*a) + T 0i ftd k7ij + T?d k p + \^d kll3 , 

Above g = detg a b = —a 2 j with 7 = det'yy. From straightforward calculations the stress-energy tensor is 



(24) 



(25) 



(26) 



T 



on 



phW 2 - p 



T 

rpO 



phW 2 (v l 



+ n 2 



(27) 
(28) 



Ri fjj ft* R3 

p hW 2 (v l - B-)(v> - !-)+p(Y 3 - ^-§-)(?9) 



phW 2 



(30) 



Note that both the fluxes and the source terms depend 
also on the primitive variables w = {p, p,e,v 1 }, and 
the source terms do not depend on derivatives of T ab . 
Eq. ([2"U)l above conserves exactly the rest-mass (or bary- 
onic mass), 



M = / d 3 x q° = / d 3 x y^/D 



(31) 



Standard York- ADM matter variables are easily recov- 
ered 



by 
Ac 
A, 



av x - p x 
n 



c s y (1 — v 2 )[y xx (l — v 2 c 2 ) — v x v x (l — c 2 )] 



(35) 
(36) 

f3 x . 



The others are obtained by permutation of indexes. 
Above the sound speed is defined by 



X + 7 K 



\ 



dP 

~dp~ 



1 

h ' 

dP 

a7 ' 



(37) 
(38) 



Note the notation conflict between x an d the metric vari- 
able defined in Scc. lII Al which is resolved in the following 
by the context of the discussion. In the Sec. lIIII we discuss 
the modeling of NS interior, i.e. the EoS of the fluid. 



III. DESCRIPTION OF NS MATTER 



PADM 

S'adm 

"ADM 



n a n b T ab = phW 2 ~p = r + D , (32) 
-n a f b T ab = phW 2 v l = S i , (33) 
l la l 3h T ah = phW 2 v l v3 + f' J p . (34) 



The system in Eq. ([2U]) is strongly hyperbolic provided 
that the EoS is causal (the sound speed is less than the 
speed of light) [65| . Eigenvalues (in direction x) are given 



The exact nature of the internal structure of a NS is 
unknown. The standard picture assumes that the matter 
of an isolated NS in hydrostatic equilibrium is strongly 
degenerate and at thermodynamic equilibrium. Conse- 
quently temperature effects are neglected and the mat- 
ter is in its ground state: cold catalyzed matter. Un- 
der these conditions the EoS has one-parameter charac- 
ter [74j |. p — P(p)- Note that if one has an one-parameter 
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EoS the GRHD equation for r is equivalent to that for D 
and thus redundant. A simple cold EoS often employed 
in NR simulations is the polytropic EoS 



P(p) = Kp T , 



(39) 



where K and T are parameters. The polytropic EoS de- 
scribes isentropic fluids and it is equivalent to the well 
known ideal gas 



P(p,e) = (T-l)pe, 



(40) 



if the flow remains smooth, i.e. no shock-heating. The 
parameter for a relativistic Fermi gas is T — 4/3, and 
1 < r < 3 for NS simulations. Note that it coincides 
with the adiabatic index of the fluid, 



f 



(41) 



The constant K in Eq. (|3"§]) fixes the entropy. Ideal gas 
and polytropic EoSs provide a simple and analytical de- 
scription of NS matter, even though it is quite rough and 
approximate. 

The NS structure consists qualitatively [75[ of an outer 
region (outer crust) that extends until the neutron drip 
Odrip ~ 10 11 gem -3 , an inner crust up to nuclear densi- 
ties, e nuc > 10 14 gem -3 , which characterize the central 
core. The composition of the inner and outer crust is rea- 
sonably well understood. Matter in the outer crust con- 
sists of a Coulomb lattice immersed in an electron gas; 
the pressure is dominated by electron pressure. As den- 
sity increases, the pressure contribution from neutrons 
becomes larger, and in the inner crust the EoS softens 
due to the attractive long-range behavior of the strong 
interaction. Models for the outer and inner crust are for 
example the BPS [z| and the BBP [z3, or HP94 [z3 
EoSs, respectively. At densities c > c nuc nuclei are not 
stable and a plasma of nucleons dominates the pressure. 
The EoS stiffens due to the repulsive short-range char- 
acter of strong interactions. The modeling of the matter 
in the core is difficult and requires assumptions on the 
nucleon-nucleon potential and the solution of the many- 
body problem. Further complications are the presence 
of hyperons and the necessity to solve the relativistic 
many-body problem, super-fluidity, pion condensation, 
and phase transition to quark matter. 

The NS core contains most of the mass (98%), while 
only a very small fraction is in the outer crust (10 -5 %). 
Despite this the (outer and inner) crust plays the major 
role in defining tidal the deformation-disruption (mass- 
shedding limit) during evolutionary scenarios. Several 
EoS for the core are proposed in the literature, see 
e.g. [79l.l80j. We refer to the set of these EoSs generically 
as realistic EoS. They are provided by tables, or, alterna- 
tively, they can be phenomenologically parametrized by 
piecewise polytropes [8lj|. Most of the realistic EoSs are 
considered "equivalent" since they are able to produce 
NSs with mass and radii that agree with observational 
constraints. See e.g. [83 - l84| and references therein for a 



detailed report on observational constraints on NS EoSs. 
Gravitational waves emission from pulsating NSs repre- 
sents a unique opportunity to constrain the EoS of the 
core 



In a dynamical scenario (e.g. merger or collapse) NS 
matter can be heated, thus acquiring a thermal compo- 
nent. This situation can not be modeled if cold (one- 
parameter) EoS are employed and the use of a hot (tem- 
perature dependent) EoS is necessary. The simplest hot 
EoS is the ideal gas Eq. (|4U|) . A simple way to extend nu- 
clear cold EoS to include a hot component is presented 
in !36|. Given values of p and e, a hot part is allowed 
evolving also the r equation and defining a "hot internal 
energy" as the difference between the actual e and the 
£ coid f rom rea li s tic cold EoS, 



£ hot = e _ e c 



(42) 



The pressure is augmented with a thermal component 
which is modeled via a simple ideal gas EoS, 



P(p,e) = P cold (p) + P hot (p, £ ) 

= P cold (p) + (r - l)pe hot 

The following relations hold: 



2P cold , d 2 e cold 



X = 



dp 2 



+(r- i)e hot - (r-i)p 
= (r-i)p. 



de 



cold 



dp 



(43) 
(44) 



(45) 



(46) 



Note that the adiabatic index in this case does not coin- 
cides with the sum of the "cold" adiabatic index and the 
ideal gas T. Hot EoSs constructed in this way are called 
hybrid. A study of the validity of this approach can be 
found in [37} . Note that obviously the hybridization of 
the polytropic EoS is the ideal gas EoS. 

The only complete (including temperature depen- 
dence) microphysical EoSs have been developed by Shen 
ct al. [Sft, Lattimer and Swesty (87}, and recently by 
Shen G. et al [88|, |89[. They allow inclusion of neutrino 
emission schemes, e.g. [9l| , and currently provide the best 
model to describe high density NS matter. 

In addition to the analytic ones, our code can handle 
cold EoSs provided by tables and implements the hy- 
bridization procedure described above. While the code 
is ready to host hot complete EoS, we do not consider 
them in this work and postpone their use to the future. 



IV. NUMERICAL METHOD 



Our code is part of the BAM code [62| - |64| . extending it 
with a module for GRHD. In the following we focus only 
on the matter solver, referring to [132] for a description of 
the code infrastructure and the algorithm for the solution 
of the Einstein equations via the BSSNOK scheme. We 
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mention only that the evolution algorithm is based on the 
method of lines (MoL) with explicit Runge-Kutta meth- 
ods (3rd order in this work) and finite differences in space 
(4th order in this work). Mesh refinement is provided by 
a hierarchy of cell-centered nested Cartesian grids and 
Berger-Oliger time stepping. Metric variables are inter- 
polated in space by means of 4th order Lagrangian poly- 
nomials and matter conservatives by a 4th order WENO 
scheme (see below for details about the latter). Inter- 
polation in Berger-Oliger time stepping is performed at 
2nd order. Some of the mesh refinement levels can be dy- 
namically moved and adapted during the time evolution 
according to the technique of "moving boxes" . GWs are 
extracted using the tp4 formalism (see Sec. Ill of 62]). 

The algorithm implemented for the matter is a robust 
high-resolution-shock-capturing (HRSC) method [9lL H3 
based on a central scheme for the numerical fluxes. It has 
been successfully used and tested in spherical symmetry 
in [93j]. HRSC schemes represent nowadays the state-of- 
the-art methods to solve GRHD equations since they can 
properly handle physical shocks, steep gradients and high 
Lorentz factors in relativistic plasma. In the following 
we describe in detail our implementation as well as the 
implementation of the EoS interface and the treatment 
of vacuum regions and spacetimc singularities. 



A. HRSC scheme for matter 

The GRHD equations are solved by means of a HRSC 
method 0, [92j which considers the semi-discrete form 
of the equations and combines the Runge-Kutta integra- 
tion with a cell-centered scheme for the RHS based on 
robust central schemes or simple Riemann solvers [94j. 
Both the time stepping and the spatial refined mesh are 
shared with the metric system. Our implementation fol- 
lows quite closely the algorithm presented in [95| (see 
also j^lMQ). The semi-discrete form of Eq. QZtfy is 



dq 



dt 



i-l,j,k 



F i+ i. j fej + other directions , 
(47) 



where h is the grid spacing and F i 



i±-k,j,k 



are the numer- 



ical fluxes (both in the x direction). For simplicity the 
description is limited to ID since all the steps necessary 
to construct the numerical fluxes can be (and actually 
are) performed in one direction at a time. In Eq. (|47[) 
the difference of the numerical fluxes is the Taylor ap- 
proximation, at a certain order r, for the divergence of 
the fluxes 

hf'ixi) = [P l+ , - 3_x) + 0(h r ) , (48) 

(r-l)/2 

F i+ i = f i+ i + £ c 2j D^f iH , (49) 

where the interface fluxes are computed with a 

Riemann solver and the D 1 * 2 ^ is a discrete operator 



which approximates derivatives of order 2j. We con- 
sider only 2nd order schemes, which amounts to drop- 
ping high-order terms in Eq. (|49j) . The interface fluxes 
are computed by the local Lax-Friedrichs (LLF) ce ntral 
scheme |94l l97j , or the well known two-speed HLL [l OQj ] 
Riemann solver. In this work we focus on LLF which, 
in spite of its simplicity, has been proved to be robust 
and co mpeti tive with respect to approximate Riemann 
also in the case of neutron star simula- 
The LLF flux reads 



solvers 101 



tions 



fi"02f 



f (LLF) 

h+i 



(50) 



The parameter a is the maximum of the local character- 
istic speeds (eigenvalues computed at interfaces i ± i) 
of the system and it is the only characteristic infor- 
mation used. The quantities written with superscript 
L/R are the physical fluxes (Eq. ([24])) and the conser- 
vative variables computed ("reconstructed") at the in- 
terface i + i. The reconstruction is performed using a 
non-oscillatory interpolation centered, respectively, on i 
{L, left) and on i + 1 (R, right). The reconstruction step 
is performed on the primitive variables. Several meth- 
ods are implemented in the code: linear Total Variation 
Diminishing (TVD) interpolation based on "minmod" 
(MM2) and "monotonized centered" (MC2) slope lim- 
iters (see e.g. |103| and jl04j ). the i nterpolat ion of the 
piecewise parabolic method (PPM) |l05l - ll07l |. and the 
third order c onvex-ess entially-non-oscillatory (CEN03) 
algorithm by [l08l Il09j . The construction of the numeri- 
cal fluxes requires an interpolation of metric quantities at 
the interfaces. We implemented both 2nd order (simple 
averages) and 4th order Lagrangian interpolation. The 
latter is the default used in all the tests presented here. 

The algorithm used to recover the primitive variables 
depends on the choice of the EoS. In case of a general EoS 
we adopt the standard algorithm described in (92|, that 
employs the EoS of the form Eq. (fl~9f and a Newton- 
Raphson method. In case of cold one-parameter EoS 
a similar procedure is adopted b ut ba sed on the equa- 
tion which defines the variable D [llOj ]. The exact equa- 
tions employed as well as the methods are detailed in 
Appendix [X] 

Boundary conditions are applied on the primitive vari- 
ables before the reconstruction step by simple extrapola- 
tion ("outflow"). 

We describe now the spatial interpolation used for the 
conservative variables in the mesh refinement, i.e. be- 
tween levels. A non-oscillatory interpolation is necessary 
in order to avoid the Gibbs phenomenon. D ifferent m eth- 
ods are adopted in different codes, see e.g. |llll ll 12j ] . As 
anticipated at the beginning of the section we ad opt a 
4th order WENO algorithm as described in [ill- The 
ID scheme is summarized in the following. Given the four 
points Xi-i, Xi, JCi+i, Xj+2 and the corresponding data 
fi—ii • fit fi+ii fi+2, two candidate interpolating polyno- 
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mials arc constructed as 



( \ £ i fi—i ( \ 
Pi{x) = fi-\ — {X-Xi) + 



fi+1 — 2fi + fi- 



-(x - Xi) 2 , 



2/i 2 

/ \ jr , — /i+2 + 4/j+l — 3/,- , . 

P2(a;) = /* H + 

/i+2 - 2/ l+ i + /j ( ,_ 2 



2fr 2 



-(a; - XiY 



(51) 



The final interpolated value is given by 

p(x) = wx{x)px{x) + w 2 (x)p 2 (x) , (52) 
where the weights are 



cti{x) 



ai(x) + a 2 (x) 

Cj(x) 
(e + IS,) 2 ■ 



(53) 
(54) 



The weights are defined in term of the smoothness indi- 
cators 

_ 25 2 64 2 13 2 

isi - + 12/i + + 

™f f --ff --f f 

TS - 25 f 2 , 64 13 2 

^2 - Y^/i + ^+1 + ^+2 + 

26 f f 52 f f 76 f t <K\ 

Y^/i+2ji — Y2 /i+2/i+1 _ Y2 ' 



and the optimal weights 

Ci(a:) 
C 2 (z) 



a;»+2 - 
3/i 

X — Xi-i 

3h 



(56) 
(57) 



In the case of a smooth function, the interpolation re- 
duces to standard 4th order Lagrangian interpolation. 
For less regular functions the order of interpolation drops 
based on the lo cal contin uity of the derivatives (see dis- 
cussions in e.g. |ll4 Ill5| ). In the implementation we set 
e = 10~ 6 to avoid division by zero. The algorithm has 
been slightly modified in order to enforce monotonicity 
in the solution [2(3] . If pi (x) at a given point is larger or 
smaller than all four function values /j, we set the cor- 
responding <Xi to zero. If all en are zero we use linear 
interpolation. 

Finally we comment about the overall convergence rate 
expected in the simulation data. The different elements 
of the algorithm described above contribute differently to 
the error budget, some errors converge away more rapidly 
while others dominate. According to our previous expe- 
rience with black hole spacetimes and to other results in 
the literature, we expect the overall error to be domi- 
nated by the truncation error of the finite difference dis- 
cretization of the right-hand-side. In this case the matter 



HRSC scheme is 2nd order accurate at most, thus the 
evolved fields and the relevant post-processed quantities, 
such as the gravitational waves and the constraints, are 
expected to converge to the continuum solution at this 
rate. 



B. Equation of state 

Polytropic and ideal gas EoSs are analytical and do 
not require special treatment. However the use of realis- 
tic EoSs, provided in form of tables, poses the problem of 
interpolation. A major requirement is that the interpola- 
tion is consistent with thermodynamics. The application 
of the first law at zero temperature translates into the 
relation 



P{p) = P 



de 
dp 



(58) 



Several methods are adopted in the literature. Widely 
used is a simple and efficient linear interpolation, see 
e.g. [116| . which is not thcrmodynamically consistent. 
Hermite polynomials can be adopted for a ther mody - 
namic ally consistent interpolation of general EoSs }ll7j |. 
After jll8l | analytic fits of the tables were also used in 
simulations (36j in a thermodynamical consistent way. 

In our code we adopt the following thermodynami- 
cally consistent procedure for cold EoS. In GRHD evo- 
lutions the quantities provided by the EoS are p, x an d 
k. They are obtained in two steps. First, for a given 
p we construct e(p) and its derivatives by interpolating 
the logarithms, i.e. the function y(x) with y = log 10 e 
and x = log 10 p. Derivatives are taken consistently from 
the interpolating polynomial. Second, the pressure is 
obtained from Eq. (|58p together with Xi which requires 
e"(p) (sec Eq. (05J). For cold EoS k = 0. Because of the 
second derivative, in order not to lose the term oc e"(p), 
a quadratic interpolation must be employed at least. We 
consider the cubic interpolation formulas given either by 
the standard Lagrangian four p oint s (centered) stencil 
or by Hermite polynomials, e.g. [l!9f . In the latter case 
a table of the derivative y'(x) must be provided, and it 
is computed consistently from Eq. (|58[). An important 
point for the accuracy of the interpolation procedure is 
the interpolation of the logarithms. In the tests per- 
formed here no relevant differences arc found between 
Lagrangian and Hermite interpolations. The method fits 
also for the cold part of the hybrid EoS described in 
SecHU 



C. Vacuum treatment 

Of particular importance in the simulation of an iso- 
lated object is the treatment of the vacuum part. The 
GRHD equations in vacuum formally do not apply and 
the numerical algorithm can not be employed directly 
because the equations to recover the primitives from 
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the conservatives are singular (see expressions in Ap- 
pendix |2| . 

The problem of simulating fluid-vacuum interfaces is 
quite generic in fluid dynamics a nd ty pically a challenge 
already at the Newtonian level |l03j . A correct, gen- 
eral and robust solution is not currently available even 
without the complication of dynamical spacetimes. A 
standard approach, largely employed in the literature, 
is instead to replace the vacuum with a minimal atmo- 
sphere of density of several orders of magnitudes smaller 
than the typical densities in the system. The main claim 
is that, since the atmosphere density is small, it will have 
a negligible dynamical impact. In case of NSs the situ- 
ation is complicated by the presence of gravity and of a 
stiff fluid. 

We implemented a simple vacuum algorithm based on 
a cold an d static atmosphere. The main ideas come 
from [23|, I llOL 1 120j ] . It consists of the following main 
prescriptions: (i) the atmosphere density value, p a tm = 
/atm max p, is chosen as a fraction, / a t m , of the maximum 
density; (ii) the atmosphere pressure and internal energy 
are chosen according to the cold part of the EoS of the 
evolved fluid; (iii) the atmosphere velocity is zero; (iv) the 
atmosphere is added to initial data in vacuum regions 
before starting the time evolution; (v) during the evo- 
lution, while recovering the primitive variables, a point 
is set to atmosphere if the density is below a threshold 
Ahr = /thrPatm- Typical values used are / atm = 1CT 10 
and /thr = 10 2 . In all our tests we found this approach 
sufficiently general and robust for our purposes. 



D. Singularity treatment 

The spacetime singularities treatment in BAM is 
based on the wel l known moving puncture method, 

e.g. [HI Bill [Ull, which relics on the BSSNOK formu- 
lation of the Einstein equations and on the gauge choice 
presented in Sec. Ill Al This method has been proved par- 
ticularly simple, elegant and robust and it is widely used 
in binary black hole simula tions flTl li"8l . lo2j as well as in 
matter simulations [27l. [47l fl22| - [l25| . It has been shown 
in fact [TIF " 



TABLE I: Initial data used in single-star evolutions. 
Columns: name, EoS, gravitational (ADM) mass M, rest- 
mass Mo, equatorial proper radius R, angular momentum J 
scaled by the square of the ADM mass, central rest-mass den- 
sity p c . Poly tropic models are computed with P = 2 and 
K = 100. 



Il27j that singularities produced by collaps- 
ing matter are naturally handled by the puncture gauge 
without particular treatment beyond standard artificial 
dissipation for the metric varia bles. 

As already pointed out in [127| however the gauge 
choice alone is not always enough to obtain stable and 
long term simulations of collapsing NSs. Unphysical val- 
ues can be produced by the HRSC scheme due to numeri- 
cal errors. They arc typically localized in a neighborhood 
of the center of the collapse and appear after the forma- 
tion of the apparent horizon. The origin of the problem 
is clear: when v — > 1 (or bigger due to numerical errors) 
the eigenvalues in Eq. ([55]) become singular and the for- 
mulas to recover the primitives, too (Appendix [A| . In 
order to prevent this our code sets the GRHD eigenval- 
ues to zero if unphysical values are computed and the 



Name 


EoS M M R 


J/M 2 


p c [xlO" 3 ] 


AO 


polytropic 1.400 1.506 9.586 





1.28 


UO 


poly tropic 1.448 1.506 5.838 





7.993 


F0 


FPS 1.400 1.566 7.367 





1.906 



numerical fluxes to zero whenever the velocity becomes 
larger than the speed of light. We found this prescrip- 
tion robust in the collapse of both a single sta r an d of 
the HMNS produced in binary simulations. In |l27| we 
tested other possibilities, in particular to set a ceiling on 
the Lorentz factor (W ce ii = 10 10 ) if the velocity becomes 
larger than the speed of light, as well as hydro-excision. 
In the case of single star collapse simulations they lead 
to comparable results, see Sec. IV Bl 



V. NS SIMULATIONS 

In this section we validate our code by considering evo- 
lutions of single NSs. Most of the tests presented here 
were suggested in [23|, Q and performed later by other 
groups. They incl ude l ong term stab ility of equilibrium 
configurations 0, llT^Jm lr28L Il29) , consistency of lin- 
ear radial oscillations 68jl, migration of an unstable con- 
figuration to a sta ble one IllOl. | l3fl . gravitational col- 
lapse to black-hole [1251 1 1281 . Il3l( 7and the evolution of 
a boosted star [23| . In the following we study the perfor- 
mance of different reconstruction procedures and the con- 
vergence of the code evolving stable equilibrium configu- 
rations. We test the EoS interpolation scheme proposed 
in Sec. IIV 51 comparing evolutions with polytropic EoS in 
analytic and table form and considering an evolution of 
a stable spherical equilibrium model described by a real- 
istic cold EoS. We demonstrate the ability of the imple- 
mented algorithm to handle shocks (migration test) and 
the formation of singularities (spherical collapse). The 
evolution of a boosted star permits the verification of a 
simple solution of Einstein equations in a fully dynam- 
ical and non-linear case and to test the moving boxes 
technique with the matter scheme. For this problem we 
additionally investigate the use of different values of the 
77 parameter in the Gamma-driver shift condition. 

The main properties of the initial data employed are 
listed in Tab. |TJ Most of them have been previously 
employed in such tests. Model AO and U0 are poly- 
tropic r = 2 spherical configurations with the same rest 
mass; AO is stable while U0 belongs to the unstable 
branch of the configurations space. Model F0 is a stable 
model computed with the FPS EoS [l3l[l33j|. Spherical 
configurations are computed with a 2-domain spectral 
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code described in |93| which employs isotropic coordi- 
nates. The values of proper radii in physical units are 
R = 14.156 km, 8.621 km and 10.879 km, respectively, 
for models AO, U0 and F0. 

The numerical grid set up for these simulations is, ex- 
cept for the migration and collapse test, such that the 
finest refinement level covers the whole star while the 
other refinement levels are used exclusively to the push 
boundary far away. This choice minimizes the effect of 
the noise due to the interpolation between boxes on the 
matter variables while it is computational more expen- 
sive for a given resolution. In the case of the migration 
and the collapse test more refinement levels are required 
to resolve the expanding or contracting fluid. Octant 
(x > 0, y > 0, z > 0) symmetry is imposed. We use the 
3rd order Runge-Kutta time stepping and a Courant- 
Friedrichs-Lewy (CFL) factor of 0.25. In some cases we 
tested also the use of 4th order Runge-Kutta finding the 
same results. All the figures in this section and in the 
following are shown for the maximum resolution unless 
explicitly stated otherwise. 



A. Stable stars 

Model AO and F0 are stable equilibrium configurations, 
therefore their evolution is trivial. Since the continuum 
solution can not evolve but has to remain in the initial 
condition, the dynamics of the numerical solution is gov- 
erned by truncation errors and by spurious effects due 
to the artificial atmosphere. This gives the possibility to 
study long term stability and convergence of the code. 
Different reconstruction schemes are considered. The 
grid is composed of three levels labeled I = 0, ..,2; the 
finest box covers entirely the star. Simulations are per- 
formed employing resolutions of h% = 0.2, 0.3, 0.4, 0.5 
(h 2 = 0.295, 0.443, 0.591, 0.738 km) for the finest box 
and last about 10 ms. 

Long term stability and convergence. We discuss here 
the simulations of model AO. It is well known, e.g. [23j ] . 
that numerical errors trigger small amplitude pulsations 
of the star which oscillates at proper mode frequencies. 
The phenomenon is depicted in Fig. Q] where the evo- 
lution of the central rest mass density normalized to its 
initial value is shown for different reconstruction schemes. 
The figure demonstrates also the ability of the code to 
maintain the initial configuration. For instance the pul- 
sations amplitude is less than 0.5% over 10 ms. Two 
frequencies dominates the pulsations: vf = 1421 Hz and 
va = 3959 Hz. They agree within the errors estimated 
from the output time sampling (2%) with the fundamen- 
tal radial linear mode and its fi r st ov ertone as computed 
by perturbative methods [134L Il35| . The figure high- 
lights a secular drift in the case that MC2 reconstruction 
is adopted. A similar drift is observed also in simula- 
tions with PPM and CENO at lower resolutions. This 
secular drift is a feature related to the evolution of the 
geometry together with the fluid since, if we perform sim- 




FIG. 1: Evolution of the central rest-mass density of model 
AO. The picture shows the normalized value p c (t)/p c (0) for 
different reconstruction schemes. 



ulations in the Cowling approximation (metric variables 
not evolved), it is almost absent at all the resolution. The 
amplitude of the pulsations is also larger in the case of 
MC2 indicating that the truncation errors are bigger as 
expected for a linear reconstruction. Notice that in the 
simulations with CENO reconstruction the overtone of 
the radial mode at frequency is more clearly visible. 

During preliminary tests with CENO reconstruction 
we observed a loss of stability between 2 and 6 ms de- 
pending on the resolution employed. We found that for 
a compact star evolution and a standard implementa- 
tion of the CENO reconstruction, the limiter tends to 
select in some points the lower order sub-stencils. A sim- 
ilar effect is also discussed in [T^|. The problem is eas- 
ily fixed in our set up by choosing a different weight in 
the weighted differences between the linear reconstruc- 
tion and the quadratic polynomial with centered sten- 
cil. Specific ally we set the free parameter a = 0.7 in 
Eq. A. 6 of [Tog to a = 0.1. In this way the limiter 
selects the central higher order stencil. We found this 
solution sufficient for all the problems and at all the res- 
olutions considered in this work. We note that also the 
PPM reconstruction has several free parameters to tune. 
We did not attem pt tu ning but instead we use the pre- 
scription given in [136]. 

Next we discuss convergence. Fig. [2] reports the evo- 
lution of the LI distance between the rest-mass density 
evolved and the initial data. The different curves refer to 
different resolutions; at every time the difference between 
two curves behaves as \\p(t) — p(0)||i oc h r + 0(h r+1 ) 
since the initial data represent also the solution for the 
evolution problem. The three panels from top left to 
bottom left refer to different reconstructions, fn all the 
cases we observe convergence with increasing resolution; 
2nd order convergence is found at early times t ~ 2 ms 
for all the methods. At late times however the MC2 
and PPM reconstructions show larger truncation errors 
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FIG. 2: Convergence in LI norm of evolutions of model AO. 
The first three panels from the top-left to the bottom-right 
show the evolution of the LI distance ||p(t) — p(0)||i. Different 
lines refer to different resolutions. Different panels refer to the 
three different reconstructions considered: MC2, PPM and 
CENO. The last panel (bottom-right) shows in a log- log plot 
the LI distances as a function of 1/h at t = 8 ms for the three 
reconstructions MC2 (solid-red line), PPM (dashed-blue line) 
and CENO (solid-green line). 



and the curves present a quadratic behavior. This leads 
to apparent over-convergent results, which indicates the 
simulations are not yet in the convergent regime. By 
contrast the use of CENO reconstruction gives 2nd or- 
der convergence over the whole simulated time. The last 
panel (bottom right) summarizes the observed conver- 
gence rate by showing in a log- log plot the LI distances 
as a function of I//12 at f = 8 ms and for different recon- 
structions. The slope of the lines gives the convergence 
rate. The over-convergent behavior for MC2 and PPM as 
well as the 2nd order convergence of CENO arc evident. 

To interpret these results we recall here that in a HRSC 
scheme the truncation errors strongly depend on the de- 
gree of smoothness of the solution and on the specific lim- 
iter employed. While the formal convergence rate of the 
methods employing different reconstruction is the same 
(2nd order), the truncation errors for this problem are 
simply different in the three cases. PPM and MC2 re- 
constructions are only first order at smooth extrema thus 
they are expected to be less accurate on smooth solution 
than CENO, but more robust in the presence of strong 
shocks. More importantly the apparent over-convergence 
behavior is not simply related to the HRSC method em- 
ployed to solve the GRHD equations, but it seems a gen- 
uine aspect of evolving the GRHD equations coupled to 
the Einstein equations. To show this we consider Fig. [31 
which is the same as Fig.[5]but in the Cowling approxima- 
tion. From an inspection of the figure it is clear that, once 
the metric is not evolved, all the reconstructions perform 
quite similarly showing perfect 2nd order convergence in 
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FIG. 3: Same as in Fig. [5] but results are computed in the 
Cowling approximation. 



LI norm at all times. Specifically one finds that CENO 
and PPM have similar performance and that for all the 
reconstructions the convergent regime is reached at the 
considered resolutions. We conclude that the slower con- 
vergence observed for MC2 and PPM in the case the 
metric is evolved is due to a combination of numerical 
errors from various part of the algorithm rather than to 
an effect related only to the reconstruction methods. 

We finally comment on rest-mass conservation. In all 
the simulations reported here we observe that the devi- 
ation runs from a maximum of AMq/Mq ~ 10 -2 with 
MC2 reconstruction at the lowest resolution to a mini- 
mum of AM0/M0 ~ 10~ 6 for CENO at the highest res- 
olution. This remarkable result is a well known conse- 
quence of the use flux-conservative methods. 

Constraint violation. The Einstein constraints are 
violated at the level of the numerical error indepen- 
dently of the formulation or the method adopted. When 
frcc-cvolution schemes are adopted, the constraints are 
only monitored (not solved) and typically the violation 
(i) grows in time, (ii) converges with increasing resolu- 
tion. In our simulations we observe this behavior. The 
biggest violation on the grid are found in the region cov- 
ered by the matter and, at least at the initial time, at 
the boundary. 

Fig. [4] summarizes our findings. The three panels from 
top-left to bottom-left show the evolution of the L2 norm 
of the Hamiltonian constraint for several resolutions in 
the cases when MC2, PPM and CENO reconstruction 
are employed. The violation converges to zero in all the 
cases. From the figures is clear that the absolute value of 
the violation is larger for the MC2 reconstruction while 
it is about a factor 10 smaller for CENO reconstruction 
when compared to MC2 and PPM. A summation of ef- 
fects and errors as described in the previous paragraph 
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FIG. 4: Hamiltonian violation during the evolution of model 
AO. The first three panels from the top-left to the bottom- 
right show the time evolution of the 2-norm of the Hamilto- 
nian constraint computed on the finest grid level which fully 
contain the star. Different lines refer to different resolutions. 
Different panels refer to the three different reconstructions 
considered: MC2, PPM and CENO. The last panel shows 
the profile in the x direction of the Hamiltonian constraint 
at time t ~ 10 ms for the maximum resolution and three re- 
constructions MC2 (solid-red line), PPM (dashed-blue line) 
and CENO (solid-green line). Note the different scales of the 
plots. 



contributes to this behavior. It is difficult to clearly iden- 
tify them due to the complexity of the equations and of 
the numerical method employed. 

We observe that, if mesh refinement is not employed 
but only one box corresponding to the 3rd level is used, 
then the norm of the Hamiltonian constraint shows an 
anti-convergence behavior at early times. This feature is 
due to the initial non-convergent constraint violation at 
the b ound ary, related to the use of Sommerfeld condi- 
tions [l37| , which dominates the violation in the interior 
at early times. The use of mesh refinement is thus impor- 
tant to minimize the effect of the boundary conditions. 

The bottom-right panel displays the spatial profile in 
the x direction of the Hamiltonian constraint at late time, 
t ~ 10 ms, for simulations employing the highest resolu- 
tion and different reconstruction schemes. The Hamilto- 
nian violation accumulates in time in the region of the 
matter. Once again the difference in the results is clear 
between the CENO reconstruction and PPM and MC2. 
In case of the unigrid simulations mentioned above the 
bulk violation dominates over the boundary violation af- 
ter the first ms of evolution in the case that PPM or MC2 
reconstruction is employed. 

These tests convinced us to use the CENO reconstruc- 
tion (see also Sec. IVI Bp . We stress that if a different 
setting is used, e.g. a different mesh, method to solve 
Einstein equations, numerical flux, etc. or simply much 



higher resolutions, the picture may change. 

The constraint accumulation and boundary condition 
effects discussed here are related to the use of BSSNOK 
and the Sommerfeld bound ary c ondition. They have 
been studied in detail in [H, Il37l | where the Z4 c for mu- 
lation was proposed as an alternative. Ref. [H, Il37l | are 
restricted to spherical symmetry, but we will investigate 
in the future the use of Z4c in 3D. 

Interpolation of EoS tables. In order to establish the 
performance of the interpolation scheme for EoS ID ta- 
bles described in Sec. IIVB1 we consider simulations of 
model AO and of model F0 employing tables. 

The evolutions of model AO performed both with the 
analytic expression for the polytropic EoS and with a 
test-table do not significantly differ. As an example the 
central rest mass density show differences smaller than 
0.01 % over 10 ms of simulated time for a resolution of 
fi2 = 0.2 (h-2 = 0.295 km). Performance is clearly affected 
by the use of tables: we observe on average a slow down 
of about 10 % but always less than 20%. The use of 
tables with different entries, e.g. 126 and 1024, does not 
lead to differences in these simulations. 

In order to obtain an accurate evolution of model F0 
we employ a t able generated from the analytic fits of the 
FPS EoS [MEIl- Initial tests with publicly available 
tables showed a progressive accumulation of round off 
errors as simulation time advances that eventually led to 
a failure during the recovery of primitives. The reason 
is the low accuracy of the table entries (in some cases 
few digits). Once tables from the fits are employed, the 
simulations are stable and accurate as in the polytropic 
case. For what concerns the proper radial frequencies of 
F0, we estimate v-p = 3045 Hz and v-r = 7232 Hz, with 
a resolution of hi = 0.2 and a total simulated time of 
7 ms. They agree wi thin few percent with perturbative 
results computed in 135|. As expected in this case the 
evolutions with the cold table and with the corresponding 
hybrid EoS do not differ significantly, while we observe 
some small differences due to the different atmosphere 
treatment. 



B. Unstable stars 

We discuss the evolution of model U0. Since it is an 
unstable configuration a small perturbation can cause ei- 
ther the migration towards a stable model of the same 
rest mass or the collapse to a black-hole. In practice trun- 
cation errors are enough to trigger the migration while 
the further addition of a small radial perturbation (bigger 
than the truncation errors) leads to the collapse. In the 
following we will consider both cases. The simulations 
presented here were performed with CENO reconstruc- 
tion. A direct quantitative comparison of the results for 
the collapse with an independent ID spherical code can 
be found in |l27l ] . 

Migration. The dynami cs of a migrating star is de- 
scribed in detail in I24l Il30l. Within the first 0.5 ms the 
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FIG. 5: Evolution of the central rest-mass density in the 
migration test. The picure shows the normalized value 
p c (t)/p c (0) in the case of an ideal gas EoS (red solid line) 
and polytropic EoS (blue dashed line). 
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FIG. 7: Dynamics of the collapse to a black hole with punc- 
ture gauge. The picture shows the evolution of the central 
rest-mass density (top panel) and of the central lapse (bot- 
tom panel). Both are normalized by the central value of the 
initial data. 



0.25 - \ 




5 10 15 20 25 30 

x(km) 



FIG. 6: Shock formation in the migration test. The picture 
shows the profiles of the specific internal energy e and of the 
component v y of the velocity at time t ~ 1.1 ms, i.e. soon after 
the shock formation. Note that only 1/10 of the numerical 
grid is shown. 



star rapidly expands while the central density decreases 
by a factor 8 to p c ~ 1.4 x 10~ 3 ; a phase of strongly 
nonlinear pulsations around the stable configuration then 
starts. The configuration reached is in fact the perturbed 
model AO, where the difference in the ADM mass be- 
tween AO and U0 has been converted into the kinetic en- 
ergy of the pulsations. If the polytropic EoS is adopted, 
thus enforcing an isentropic evolution, oscillations can 
be damped by the effective numerical viscosity of the 
HRSC scheme and by spurious interaction with the at- 
mosphere. If the ideal gas EoS is employed, thus allow- 
ing shock heating, a shock forms during the first pul- 



sation at the interface between the inner core and the 
infalling mantle (outer lower density matter). The in- 
ner core then bounces and expands again several times 
feeding the shock with kinetic energy, which dissipates 
it into thermal energy. As a result the oscillation ampli- 
tudes decrease. Fig. [5] shows this behavior captured by 
our simulation; the evolution of central rest-mass density 
is plotted for the two different EoS. The amplitude of the 
pulses remains approximately constant during the simu- 
lation in the case of polytropic EoS, indicating the spu- 
rious numerical effects arc negligible, while it decreases 
as expected in the case of the ideal gas EoS. The capture 
of the shock dynamics is demonstrated in Fig. [5] which 
shows the profile in the x direction of the specific internal 
energy and of the x component of the velocity at the time 
£ ~ 1.1 ms, i.e. just after of the shock formation when 
the shock wave is moving outward. The figure clearly 
indicates how part of the matter is expanding outwards 
while the low density mantle is infalling towards the sym- 
metry center. Note the steep profile of the velocity. We 
observe that, in the ideal gas EoS case, the expansion of 
the star extends to the boundary of the numerical grid 
and a fraction of the mass, 1-2 % percent over the simu- 
lated time, is lost outside |168| . The grid covers a cube 
of side 200 (295 km) and five levels are employed with a 
maximum resolution of hi = 0.15 (hi = 0.221 km) and 
a minimum of ho = 2.4 (ho = 3.544 km). The region 
within a coordinate radius of about r ~ 50 (74 km) is 
entirely resolved by the two finest refinement levels. By 
contrast, in the isentropic case, the expansion reaches a 
maximum coordinate radius of r ~ 25 (37 km) and the 
rest mass is conserved as in the stable star tests. 

Collapse. The collapse is triggered by introducing a 
radial perturbation of the velocity field with an ampli- 
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tude larger than truncation errors. Since we are only 
interested in testing the ability of the code to handle the 
formation of singularities, to keep the set up simple we 
do not solve the constraints after imposing the perturba- 
tion. In our experience this procedure does not introduce 
relevant unphysical effects, while it is clearly an inconsis- 
tent way to solv e Einstein equations (see also discussions 
in (l9llllCtll38{ ). Furthermore the re sults obtained are in 
close qualitative agreement with [24], Il30j where the con- 
straints were solved. For these simulations we employ 
an eight level grid with maximum resolutions of — 
0.05, 0.03125, 0.025 (h 7 = 0.0738, 0.0461, 0.0369 km) 
which cover a cube of side 1.164 (1.72 km). 

The collapse happens in the first 0.3 ms of simula- 
tion: the matter falls towards the symmetry center while 
the metric varies rapidly adapting itself to the mat- 
ter distribution. At about £ah ~ 0.23 ms an appar- 
ent horizon forms with an initial coordinate radius of 
r AH («AH) < 0.96 (1.41 km or 0.66 M). Part of the mat- 
ter is outside of it and then rapidly accreted. Fig. [7] 
shows the evolution of the central rest-mass density (top 
panel) and that of the central lapse (bottom panel). The 
central density increases of a factor two at £ah while 
the lapse collapses towards zero. After £ah the gauge 
conditions in Eq. (fT4")) and Eq. (fT5|) . in particular the 
Gamma-driver for the shift, play the main rol e in han- 
dling the singularity. As described in detail in |l27| . the 
shift condition deforms the spatial coordinates pushing 
the proper (Schwarzschild) radii rsd™ ^ 1-7 (2.55 km 
or 1.2 M) outside the innermost point of t he n umerical 
grid (sec also the top panel of Fig. 6 in }l27l ]h Con- 
sequently: (i) matter effectively disappears from the nu- 
merical grid; (ii) the end state of the collaps e is the trum- 
pet slice of Schwarzschild [ll, [HI EH, . The central 
rest-mass density in Fig. [7] reaches the atmosphere value 
about At ~ 0.1 ms after £ah- Some spurious effects are 
visible during this time interval and they are related to 
the numerical treatment of the matter fields discussed in 
Sec. IIV Dl We experimented with the alternative treat- 
ments described in Sec. IIV D[ finding unimportant quan- 
titative differences only in the short phase before all the 
matter disappears except for the atmosphere. 

The behavior of the matter is also demonstrated in 
Fig. [5] which shows the evolution of the irreducible mass 
of the black hole normalized to the ADM mass of the 
system (solid line) and of the rest mass normalized to 
the initial value (dashed lines). The irreducible mass of 
the black hole is obviously zero at the beginning and 
after £ah it rapidly reaches a value corresponding to the 
initial ADM mass of the system. The rest mass is shown 
computed on the finest refinement level (level 7) for the 
best resolution (hi = 0.025) and on a coarser one (level 5) 
with a resolution of h^/A. Since the matter is contracting, 
the mass on the level 7 initially increases and after £ah 
the level contains all the matter. Level 5 instead always 
contains all the matter but after £ah does not resolve it 
properly. During the collapse the rest-mass is conserved 
up to 0.05% while, after £ah, h rapidly decreases to the 
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FIG. 8: Evolution of different mass quantities in the collapse 
test. The picture shows the irreducible mass (solid line) of 
the final black hole and the rest mass of the matter (dashed 
line) normalized to their initial values at level 5 and 7. 



atmosphere value, M (t > t A n) ~ 10~ 6 il/ (0). 



C. Boosted star 

We discuss here the evolution of model AO boosted 
in the x direction at a speed of v = 0.5 corresponding 
to a Lorentz factor of W = 1.16. The setup of the 
initial data is as described in [23|. The test is inter- 
esting because it gives the possibility to experiment in 
a simple scenario with several points: fully dynamical 
and nonlinear evolutions, the performance of the moving 
boxes with matter, and different gauge conditions. We 
use five refinement levels, the moving one is the finest 
which entirely covers the star. The different resolutions 
employed for the finest box are hi = 0.4, 0.278, 0.208 
(h 4 = 0.588, 0.408, 0.306 km). 

The solution of the evolution problem depends on the 
gauge conditions employed. If we choose to simply advect 
lapse and shift, the solution is analytic and it is just a 
time shift of the initial metric and matter profiles. If the 
conditions in Eq. (fT4")) and Eq. (fl"5|) are used, the solution 
is not analytic. In order to investigate the numerical 
solution obtained under the 1+log and Gamma-driver 
condition we consider evolutions with different values of 
the parameter r\. A similar investigation has been carried 
out for BNS and it is presented in Sec. IVIBI Fig. [5] 
summarises our findings. It shows the profiles of the rest- 
mass density in the x direction at the final time of t = 
1.72 ms for evolutions using different values of rj. All the 
profiles are shifted back to the initial position of the star, 
the initial profile is also plotted. As apparent from the 
figure, the choice 77 = corresponds to the case closest 
to the analytic solution, while for higher values the star 
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FIG. 9: Profile of rest-mass density of the evolved boosted 
model AO. The picture shows the rest-mass density profile 
normalized to the initial value at time t = 1.72 ms for evolu- 
tions corresponding to different values of the 77 parameter in 
the Gamma-driver shift (dashed colored lines). The profiles 
are shifted back to the initial position of the star, the initial 
data is also plotted (black solid line). 
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FIG. 10: Self-convergence test of the evolution of a boosted 
star. The pictures shows the differences between the pro- 
files of the rest-mass density at t — 3.95 ms evolved with 
different resolutions. The labels of the solid lines "low", 
"med" and "hig" in the legend correspond to resolutions 
hi — 0.4, 0.201, 0.278. The difference between medium and 
high resolution is scaled for 2nd order (dashed line) and 1st or- 
der (dotted line). The gamma-driver condition employs rj = 0. 



profile is progressively more distorted in the direction of 
motion. 

We finally comment about convergence. Fig. [TU] shows 
the point-wise self-convergence the spatial profile of p at 
t = 3.95 ms for 77 = 0. We show the differences between 
the rest-mass density profile computed at different reso- 
lutions. The difference between the medium, /14 = 0.278, 
and the high resolution, hi = 0.208, is scaled by a factor 
corresponding to 1st (dotted line) and 2nd (dashed line) 
order convergence. Point- wise convergence is lost at early 
times but the magnitude of the errors scales at about 2nd 
order thus indicating convergence in the LI norm as ex- 
pected. Note that the evolution time presented here is a 
factor of about 10 4 longer than simulations in [23| . 



VI. BNS SIMULATIONS 

In this section we discuss our BNS simulations. We 
focus on a simple equal-mass irrotational configuration 
evolved with polytropic and ideal gas EoS. In both evo- 
lutions the formation of a HMNS is observed, but in the 
isentropic case it collapses after about 3 ms producing a 
Kerr BH surrounded by a disk of mass Md < 2 x 10 _2 Mo, 
while in the other case the HMNS survives for about 9 ms 
before collapsing. We give an overview of the dynam- 
ics and discuss the impact of using different resolutions 
and different reconstructions schemes. We present results 
concerning the use of different values of the damping pa- 
rameter 77 in the shift condition: while small values of 
77 lead to less coordinate eccentricity during the inspi- 
ral, they reduce the coordinate size of the final BH. Wc 



TABLE II: Parameters of the initial binary configuration. 
Columns: name, ADM mass M, rest mass Mi and ADM mass 
of each star in isolation, angular momentum J scaled by 
the square of the ADM mass, gravitational wave frequency wo, 
proper separation d, central density of each star p c . The pa- 
rameters in the polytropic EoS are r = 2 and K — 123.6489. 

name M M h M+ J/M 2 ujq [Hz] d p c [xlO~ 4 ] 

G2P14 2.998 1.625 1.499 0.4450 589 36.582 9.569 



present and characterize the GWs computed from the 
simulations. We compute the actual GWs degree of free- 
dom (metric waveforms) by means of the post-processing 
algorithms described in [671 |68| and [69j and compare 
their performance. Convergence tests are performed in- 
dicating 2nd order convergence of the inspiral waveforms 
without any time-shifting procedure. For the first time in 
BNS simulation, we estimated precise error-bars on the 
waveforms by extrapolating the results in resolution. 

A. Initial data and numerical settings 

We employ as initial data quasi-equilibrium configura- 
tions of irrotational equal-mass binaries in quasi-circular 
orbits. These initial configurations are computed with 
a multi-domain spectral code which solves the Einstein 
constraint equations under the assumption of a confor- 
mally fl at m etric. The code is based on the LORENE 
library jl4ll | and provided by the NR group in LUTH 
(Mcudon). These initial data represent to date the most 



15 



accurate computation of equilibrium BNSs and they are 
publicly available on the web. The assumption of irro- 
tational flow is believed to be astrophysically realistic 
since the spin period of the neutron star is larger then 
the orbital period in the late inspiral. In other terms, 
the time of coalescence due to gravitational radiation is 
shorter than that of synchronization due to viscosity. Be- 
cause the late inspiral and merger phases arc expected to 
happen a very long time after the birth of the binary 
system (~ 10 s yr), the temperature of the matter can 
be assumed to be negligible. The neutron stars in the 
equilibrium system are thus described by the cold EoS. 

The properties of the initial configuration that we con- 
sider in the following are summarized Tab. [TTJ The bi- 
nary has ADM mass M = 2.998, angular momentum 
J = 8.836 and proper separation d ~ 36.582 (54 km), 
thus the compactness of the system is M/d ~ 0.08. The 
coordinate separation is 30.457 (45 km). The rest-mass 
and ADM mass of each star in isolation (d — > oo, spher- 
ical configuration with same rest-mass) are Mb = 1.625 
and M* = 1.499, respectively. Note the notation for the 
rest-mass of the star in isolation, Mb, and for the rest- 
mass of the binary, Ma. The compactness of each star in 
isolation is M+/R = 0.14. The equilibrium configuration 
was first computed in [66|, Il42l | . 

As a test for the implementation we performed evolu- 
tions with other initial configurations. In particular we 
considered an irrotational equal- mass bina ry described 
by the APR EoS [H3| computed in [HI, [HI and evolved 
in [2f| Il45l | . The evolutions were performe d usi ng both 
an APR table based on the analytic fit of [l45l | and its 
hybridization with a T = 2 ideal gas EoS. The robust- 
ness of the procedure described in Sec. IIVBI has been 
checked also in the case of BNSs. All these tests were 
performed at grid configuration L, M and HO (see be- 
low), thus we do not discuss them in the following. An- 
other configuration whose evolution was explored is an 
irrotational equal-mass binary described by a polytropic 
EoS and with ADM mass M = 3.2 51 and compactness 
of each star = 0.16 0, [HI- As expected evo- 

lutions with both polytropic and ideal gas EoSs led to 
a prompt collapse to a BH and the main evolution hap- 
pens in 11 ms. Results are comparable to the findings 
in [6(|. Since the evolution of model G2P14 has a more 
complicate post-merger dynamics we focus on the results 
obtained from these evolutions. They evidently represent 
a more stringent test for the code. 

In Tab. IIHI we report the grid configurations used for 
the evolution simulations. Simulations with configura- 
tions L and M can be run on a small machines. They 
need between 8 and 16 processors with 1GB of memory 
per core. While they can be carried out without any 
problem, thus proving the flexibility of the code, the re- 
sults are too inaccurate to be considered for a sensible 
analysis. Convergence can be measured only in simula- 
tions employing configurations HO as shown in the fol- 
lowing. For each grid configuration the finest refinement 
level covers each star entirely. The latter is resolved with 



TABLE III: Grid configurations used in BNS simulations. 
Columns: identifier, number of grid levels (boxes), number 
of moving boxes, resolution of finest box (dimensionless and 
km), number of points in finest box, resolution of coarsest box 
(dimensionless and km), number of points in coarsest box. 

name RL MRL hs h$ [km] N$ ho ho [km] No boundary 



L 


6 


4 


0.500 


0.74 


40 


16.0 


23.6 


80 


945 


M 


6 


4 


0.400 


0.59 


50 


12.8 


18.9 


100 


945 


HO 


6 


4 


0.313 


0.46 


64 


10.0 


14.8 


128 


945 


HI 


6 


4 


0.250 


0.37 


80 


8.0 


11.8 


120 


709 


H2 


6 


4 


0.200 


0.29 


100 


6.4 


9.4 


210 


982 


H3 


6 


4 


0.156 


0.23 


128 


5.0 


7.4 


260 


960 



128 points in runs H3 and 64 points in HO. Let us briefly 
comment about the grid setti ng used h ere and those used 
for BBH simulations, e.g. [58L l62l fl46l ] . For BBH simula- 
tions roughly half the number of points per direction than 
here is used together with more grid levels (typically 9-11 
levels). Higher resolutions are reached in the finest grid 
level in order to resolve the punctures, while comparable 
(or less) resolutions are used on level 5. Therefore the 
horizon of the final BH in BBH simulations is resolved 
typically on level 6 or 7 with a resolution about two to 
four times better than in BNS simulations. One impor- 
tant consequence is that the precision of the apparent 
horizon finder is affected. 

The performance of the code for each grid configura- 
tion is reported in Tab. IIVI The BNS runs described 
in this paper require an average speed of ~ 3 M/hr 
(H3) on 128 processors in the LRZ cluster (Munich) and 
~ 9 M/hr (HO) on 32 processors. In physical units 10 M 
of the configuration selected corresponds to ~ 0.05 ms of 
simulation. We performed scaling tests up to 512 pro- 
cessors finding good scaling properties by using larger 
grid setups and higher resolutions. Production runs with 
resolutions of hs ~ 0.12 — 0.10 on 256-512 processors are 
thus definitely feasible in reasonable times with our code. 
Here we did not run such kind of simulations because of 
our computer time restrictions. 

Our grid settings are similar to those of other codes (22l . 

Ill2j . The highest resolution employed here is 30 % 
lower than the maximum resolution used to date on 
BNS simulations employing mcsh-refined-Cartesian-grid- 
based codes [22|, [56| and comparable to [2(| . 

If not explicitly stated the data discussed refer to simu- 
lations employing RK3 time stepping, CFL factor of 0.25 
and CENO reconstruction. 



B. Overview of the dynamics 

The initial configuration has been evolved using both 
the cold polytropic EoS and the ideal gas r = 2 EoS. In 
the following we will refer to the isentropic evolution with 
the name of the initial configuration while we will use the 
suffix "hot" when thermal effects arc simulated. Evo- 




FIG. 11: Dynamics of the G2P14 evolution. The picture shows contour plots in the x — y plane of the rest-mass density p and 
the velocity v* at different times. Data refer to run H3 
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TABLE IV: Performances of BNS runs. Columns: grid con- 
figuration, number of processor, cpu memory, total runtime, 
average speed. Runs last to t — 5000(1666 M) on LRZ Mu- 
nich. The numbers include LORENE initial data interpola- 
tion and checkpointing. 



Grid 


nproc 


mem [Gb] time [CPU hr] 


speed [M/hr] 


HO 


32 


90 


192 


8.7 


HI 


32 


96 


268 


6.2 


H2 


128 


120 


254 


6.5 


H3 


128 


165 


480 


3.5 



lutions of this configurations were previously performed 
in [22l. |27l| . Overall, we observe a qualitative agreement 
between our simulations and the results reported in the 
literature. Main quantitative differences are found in the 
post-merger phase, in particular in the collapse time of 
the HMNS. This is not very surprising since soon after 
the two stars come in contact the convergence of this 
kind of simulations drops to first order, e.g. |(3(|. The 
post-merger phase is thus very dependent on the reso- 
lution and grid settings employed as well as on the spe- 
cific HRSC scheme employed. In the following we will 
mostly concentrate our discussion on G2P14 evolutions 
since they were performed at all the grid configurations. 
The G2P14hot evolutions were performed up to HI con- 
figuration. Comparison between G2P14hot and G2P14 
are presented for HI. 

Fig. Qj] shows contour plots in the x — y (orbital) 
plane of the rest-mass density and the velocity field at 
different times of the G2P14 evolution. Data refer to 
run H3. The binary performs about 3 orbits before 
the merger. We define the merger time t m as the peak 
of the (2,2) mode of the GW amplitude, I/122I, where 
h = h-\- — 1 h x (see next section). Clearly the two stars 
come in contact before the merger time. From the H3 
run we have t m — 1765 (8.69 ms) while the contact time 
is about 1290 (6.3 ms). After the merger we observe 
a bar-shaped differen tially rotating star with rest-mass 
~ 2M b : the HMNS pJ69j |. Note in Fig. HU the initial 
rotational symmetry of the HMNS (obtained without 
imposing 7r-symmctry in the grid) and how the sym- 
metry is broken during the following dynamics. The 
large non-axisymmetric d eformation of the HMNS causes 
a strong GWs emission jl47l . Il48j which carries away 
matter angular momentum. As a result the HMNS 
becomes more compact and finally collapses at about 
^ah ~ 2118 (10.43 ms) when an apparent horizon is 
first formed. A fraction of the matter remains outside 
the horizon and forms an accretion disk. Note that in 
the disk rotational symmetry is approximately restored. 
The mass and spin of the BH estimated from the appar- 
ent horizon are rather inaccurate due to a lack of res- 
olution. They first rapidly grow in time reaching local 
maximum, then the BH mass is observed to increase (see 
Fig. [T9]) while the angular momentum decreases. We 
estimate at the end of the evolution Mbh ~ 2.77 and 




t(ms) x(km) 

FIG. 12: Orbital motion of G2P14. The figure shows the 
evolution of the coordinate (top-left) and proper (bottom- 
left) separation of the binary and the star-track (right) of the 
two stars for different grid settings. 



&bh = Jbh/M# h ~ 0.72. They have discrepancies of 
about 7% and 25% with the expected value once the 
gravitational radiation emission has been taken into ac- 
count. Note that the collapse simulation discussed in 
Sec. IV Bl where the apparent horizon is computed with 
satisfactory accuracy, employs in the finest grid level a 
resolution about ten times higher than in the HI runs. 
Additionally the gauge choice plays an important role in 
determining the coordinate size of the BH as discussed 
below. 

The mesh refinement implementation is such that as 
soon as the two moving boxes come in contact (which 
happens before the contact of the two stars) a larger 
box with the same resolution is constructed. Before the 
two stars come in touch however it can happen that the 
boxes split back to the initial ones, evolve individually 
and merge again. The reason is that the evolution of a 
very large box is computationally not affordable in terms 
of memory and time. This behavior is well tested in BBH 
simulations [62j . It could lead to a lack of accuracy at the 
center of the grid but in practice it has a negligible im- 
pact since it happens when the main part of the matter 
is still distributed away from the center. 

In the right panel of Fig. [T^] we report the star-track 
of one of the stars, defined as the minimum of the lapse 
for resolutions HI, H2 and H3. The left panel instead 
shows for the same resolutions the coordinate (top) and 
proper (bottom) spatial separation. As is evident from 
the figure, the orbital motion has some eccentricity due 
to the initial data (visible in the proper separation plot) 
and a coordinate eccentricity (visible in the star-tracks 
and as oscillations in the coordinate separation) due to 
the evolution itself. The eccentricity of the initial data 
is caused by the conformally-flat approximation; the ef- 
fect becomes bigger at smaller separation. In [26| the 
problem is corrected by adding an "approaching" radial 
velocity based on the post-Newtonian T4 formula. Since 
the procedure is constraint violating we prefer not to 
adopt it. The contribution to the coordinate eccentric- 
ity is mainly related to the shift condition and discussed 
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FIG. 13: Evolution of the central rest-mass density in BNS 
simulations. Runs G2P14 are reported for the three highest 
resolutions while run G2P14hot for HI. 
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FIG. 14: Hamiltonian violation in G2P14 evolution. The 
figure shows the evolution of the L2 norm of the Hamiltonian 
constraint for different resolutions. 



below. The use of a lower resolution results in an earlier 
merger. While the truncation errors due to finite resolu- 
tion are quite large, the simulations with grid settings H 
have entered the "convergent regime" in the sense that 
we are able to estimate convergence. 

The evolution of the maximum rest-mass density is 
shown in Fig. [I]J] for different resolutions. After the 
merger about two quasi-radial oscillations of the HMNS 
are observed. The maximum density increases by about 
a factor two indicating that the HMNS becomes more 
compact until, finally, it collapses. The use of more res- 
olution results in a more persistent HMNS (see below 
for a comparative discussion with G2P14hot). Note the 
quite large differences between the runs at different reso- 
lutions. After the collapse the matter of the disk accretes 
onto the BH. 

Concerning the constraints, the biggest violation is 
found in the Hamiltonian. The momentum constraint 
is generically one order of magnitude smaller and be- 
comes comparable to the Hamiltonian constraint only 
after the formation of the "puncture" . The evolution of 
the L2 norm of the Hamiltonian constraint is displayed 
in Fig. 1141 The preservation of the constraint improves 
with resolution. Wc find 2nd order convergence during 
the whole inspiral. Most of the violation is observed in 
the region covered by the matter, similarly to what was 
discussed for the test involving a single star spacctimc. 
After the two stars come in contact the violation rapidly 
increases and the convergence rate drops down to first 
order. At the time when an apparent horizon forms the 
norms show a large gradient. 

The conservation of the rest-mass is excellent until the 
BH forms with a largest deviation of AMq/Mq ~ 1 % 
in the HI runs and a minimum of AMq/Mq ~ 0.5 % 
in the H3 runs. This behavior is illustrated in Fig. Q~5] 
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FIG. 15: Conservation of rest-mass in BNS simulations. 
The figure shows the evolution of the rest-mass normalized 
to the initial value for G2P14 and G2P14hot evolutions. 
Runs G2P14 are reported for the three highest resolutions, 
G2P14hot for resolution HI. 



which displays the evolution of AMo /Mq. The plot shows 
also matter accretion after the collapse and indicates an 
upper limit (from run H3) for the initial rest-mass of the 
disk, M<j <2x 10 _2 Mo, since the integral in Eq. (j3Tj) is 
performed on the whole grid (including the BH interior). 

We finally compare the dynamics of G2P14 with 
G2P14hot. The inspiral motion of the two binaries 
present small differences due to spurious numerical ef- 
fects. Small errors, triggered by the artificial atmosphere 
treatment, propagate as simulation time advances and 
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FIG. 16: Shock formation in G2P14hot evolution. The figure 
shows the quantity K = p/p F normalized by its initial value 
on the orbital plane at time t — 1275 (6.256 ms) in a log 10 
scale. Data refer to run HI. 



artificially heat the stars especially at their surface. The 
effects of these errors on the GWs are quantified in the 
following section. When the two stars touch physical ef- 
fects dominate instead and the evolutions significantly 
differ. Fig. 1161 shows a contour plot in the x — y plane of 
the quantity K = p/p T (normalized by its initial value) 
for the G2P14hot evolution. K gives a simple measure 
of the entropy wh ose v ariation indicates the presence of 
a shock (see e.g. j 1491 ] for the development of a "shock 
detector"), and in the G2P14 evolution K is constant 
by construction. In contrast, during G2P14hot evolu- 
tion the inclusion of thermal effects yields shock forma- 
tion when the arms of the two distorted stars come in 
contact. The figure shows exactly this phenomenon indi- 
cating that the quantity K increases by about 5 orders 
of magnitude. The thermal energy of the fluid rapidly 
increases reaching peaks of e hot - 0.025 - 0.03 corre- 
sponding to temperatures of T ~ 2 x 10 11 K. At the time 
of the collapse the average temperature of the fluid is of 
order T ~ 6 x 10 10 K and it has reached peaks up to 
T~2x 10 12 K. 

Due to the additional pressure support provided by 
thermal effects (see Eq. the HMNS in the G2P14hot 
evolution collapses later. Fig. Q2] indicates that the col- 
lapse takes place after ~ 9 ms after the HMNS formation. 
During this time interval several quasi-radial oscillations 
at a frequency vp ~ 473 Hz and a strong GWs emis- 
sion are observed (see next section) . The central density 
shows again a drift to larger values indicating that the 
star is becoming progressively more compact. Results 
similar to G2P14 are found for what concerns rest-mass 
conservation as demonstrated in Fig. 1151 

In (2?| the evolution of G2P14 leads to a prompt col- 
lapse without the formation of a HMNS, while the evolu- 



tion of G2P14hot leads to a collapse at about t ~ 14 ms 
(si mula tions employ the PPM implementation described 
in [l50l | and the Marquina numerical flux). In [22j], where 
a different grid setting is employed (but the same code 
with the HLL numerical flux), the evolution of G2P14hot 
leads to a collapse at about t ~ 16 — 17 ms. In both 
works [12, H3, and differently from here, the rotational 
binary symmetry (7r-symmetry) is enforced in the nu- 
merical grid but this procedure is expected at most to 
lead to a more persistent HMNS due to non- linear mode 
couplings (in particular m — 1 modes) |l48l |. It is thus 
unlikely that the grid symmetry is the reason for the 
discrepancy as demonstrated below. In the following we 
show how the HMNS dynamics is very sensitive to the nu- 
merical method employed presenting results from simu- 
lations that employ different reconstruction methods but 
are otherwise equivalent. Together with the low conver- 
gence of the HRSC methods in presence of turbulence 27 1 
or shocks and given a dependence on the grid settings 6C | 
the differences observed are not surprising. 

Effect of reconstruction methods. As pointed out in 
Scc.[V]the different reconstruction methods have different 
truncation error magnitudes and once inserted in the al- 
gorithm produce quantitatively different results. Ref. (30T | 
already pointed out that the use of very dissipative lim- 
iters such as MM2 can lead to very different waveforms 
with respect to (formally) more accurate methods while 
maintaining the same nominal order of convergence. Our 
simulations confirm this conclusion and here we show how 
dramatic the influence of a dissipative scheme on the dy- 
namics can be. 

We consider G2P14 evolutions at resolution HI with 
the following reconstructions: MM2, MC2, PPM, CENO 
and CENO based o n th e MM2 sub-stencil (instead of 
the standard MC2 [l09j ] ) . Since the resolution is not 
optimal we expect the effects due to different trunca- 
tion errors to be strongly magnified. Fig. [T7] shows the 
evolution of the maximum rest-mass density for the dif- 
ferent methods. The vertical line marks the maximum 
of p computed from the H3 run. A visual inspection of 
the data is already conclusive in this case. As partially 
expected, the most diffusive schemes lead to an earlier 
merger and collapse while formally high-order reconstruc- 
tions are closer to our best estimate. In the simulation 
with MM2 reconstruction the inspiral is very short, the 
HMNS is not formed and the binary evolution results in 
an early prompt collapse. The simulations with PPM and 
MC2 give similar results slightly speeding up the collapse 
with respect CENO. When we introduce a more dissipa- 
tive component in the CENO limiter, i.e. the MM2 linear 
sub-stencil, the global results are significantly affected: 
an earlier prompt collapse is observed. The specific rea- 
son is obviously the fact that the linear sub-stencil MM2 
is often selected by the limiter. We mention here that, 
according to our experience, MM2 reconstruction does 
not guarantee the long-term-stability of an equilibrium 
spherical star even in ID [93|. 

Clearly the big differences shown here become progres- 
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FIG. 17: Dependence of the dynamics on the use of different 
reconstruction schemes. The figure shows the maximum of 
the rest-mass density obtained from G2P14 evolutions with 
different reconstruction methods. Reference vertical line cor- 
responding to max p of run H3. 



sively smaller when higher resolutions are considered. 
However, due to the slow (2nd) order convergence of 
HRSC, we expect they play a major role also at the higher 
resolutions employed for standard production runs. We 
point out that in comparing some runs with grid con- 
figurations L, M and HO we observed a non-monotonic 
behavior of the results for increasing resolution. This 
effect seems to be related to the poor resolution but wc 
can not exclude that the phenomenon will also happen at 
higher resolutions. Finally, the figure demonstrates that 
the use of 7T-symmetry does not have a relevant role (at 
least at this resolution) in the persistence of the HMNS 
as mentioned above. 

Gauge study: the rj parameter. Among the param- 
eters entering the gauge conditions in Eq. (fT4")) and 
Eq. (|15[) the damping parameter rj in the Gamma-driver 
shift equation has been recently the subject of in vesti - 
gati on for BBH simulations with unequal masses [151 
Il57l |. In case of equal mass BBHs it is typically set to 
r]=l/M or rj = 2/M, where M is the sum of the ADM 
masses of the two holes (M = 1), and it has been proved 
to be important to properly resolve the holes on the co- 
ordinat e gr id [62j ■ In the case of BNSs it has been sug- 



gested jl 12( | that results are not significantly affected by 
this choice. Differe nt val ues are adopted in the literature, 
e.g. 77 = 3/A4 [H EH, V = 1 Wk- without a detailed 
analysis. We demonstrate in the following that in our 
set up different choices of rj do have an influence on the 
dynamics of the inspiral. In particular smaller values of 
rj lead to a less eccentric coordinate orbital motion and 
to a smaller coordinate size of the final BH. 

We consider G2P14 simulations with grid config- 
uration HI and values rj — 0, 0.1, 0.3, 0.6, 0.9, 1.8, 
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FIG. 18: Effect on orbital dynamics of different choices of 
the parameter rj in the Gamma-driver shift condition. The 
top panel shows the star-track of one star. The central panel 
shows the evolution of the coordinate separation while the 
bottom panel shows the evolution of the proper separation of 
the binary. Each line refers to a different G2P140 evolution 
with a different rj. Data refer to HI runs. 



i.e. 77 ~ 0/M, 0.3/M, 0.9/M, 1.8/M2.7/M, 5.4/M. 
Fig. [TBI shows the star-track (top panel), the evolution of 
coordinate separation (central panel), and of the proper 
separation (bottom panel) from simulations with differ- 
ent values of rj. The coordinate separation shows large 
oscillations and a systematic shift of the merger to ear- 
lier times for higher values of rj. The latter is also visi- 
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FIG. 19: Effect on the final BH mass and coordinate radius 
of the parameter r\ in the Gamma-driver shift condition. The 
top panel shows the irreducible mass of the final BH normal- 
ized by the ADM mass of the system as computed from the 
apparent horizon finder. The bottom panel shows the coor- 
dinate radius of the apparent horizon. Each line refers to a 
different G2P140 evolution with a different r\. Data refer to 
HI runs. 
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FIG. 20: tpi waveforms from G2P14 evolution. The figure 
shows the real part (top panel) and the imaginary part (bot- 
tom panel) of the r ip22 waveform extracted at r = 200 from 
the H3 run. The amplitude is also shown in both panels. 



Gravitational waves 



ble in the evolution of the proper separation, which in- 
stead highlights non-circular effects due to the confor- 
mally flat initial data. The star-track shows how the co- 
ordinate eccentricity progressively reduces when smaller 
values of r\ are employed. When considering the gravi- 
tational wave emission the coordinate eccentricity does 
influence the extracted waves and it manifests itself in a 
phase difference accumulating during the inspiral. Tak- 
ing as a reference the r\ = case, the difference in the 
phase <j> of r 1P22 ( see Sec. IVI C|) computed with 77 = 
is of A0 = 0.43, 0.90, 2.25, 3.20 rad respectively for 
77 = 0.1, 0.3, 0.6, 0.9 at the merger time. 

Fig. [TO] shows the irreducible mass, Mah (normalized 
to the ADM mass of the system), and the coordinate ra- 
dius, tah, of the apparent horizon. As expected the co- 
ordinate size of the final BH is larger for higher values of 
77. The choice of 77 must be guided by the requirement of 
minimizing the coordinate eccentricity of the orbital mo- 
tion while keeping the ability of correctly resolving the 
final BH on the finest grid level. The use of 77 = 1.8 for 
example is not shown in the figure because it produces an 
unacceptable eccentricity in the dynamics (see Fig. I18p 
and the final BH has a size not completely contained in 
the refinement level 5, resulting in very inaccurate re- 
sults. In a similar way the use of large 77 also affects the 
coordinate radius of the GW extraction spheres; smaller 
proper radii corresponds to higher 77. 



Gravitational radiation plays the fundamental role 
driving the dynamics of the binary system. GWs en- 
code information from each phase of the evolution from 
the inspiral to the collapse and the BH-disk formation. 
We analyze in the following the GWs computed from our 
simulations. Again we focus the discussion on G2P14 
evolution. 

Our simulations show that about 1 % of the initial 
ADM energy and 13 % of the initial angular momentum 
are emitted in GWs during the simulation. The main 
emission channel is the (£,m) = (2,2) mode of the mul- 
tipolar (s = —2 spin-weighted spherical harmonics) de- 
composition of GWs. The (2, 2) mode includes about 
97 % of the entire radiated energy, thus we focus on the 
analysis of this mode. Fig. [50] shows the complex r 1P22 
waveform computed from the Newman-Pcnrose scalar 7^4 
as described in [62j from run H3. GWs are extracted 
on coordinate spheres on grid level 1. The complex 
waveform is usually decomposed in amplitude and phase, 
r " t l J 22 = ^4exp(— 1(f). Fig. I2TJ1 shows the real (blue solid 
line) and the imaginary part (red dashed line) as well as 
the amplitude (green solid line). The extraction sphere 
is placed at coordinate radius r = 200 (295.3 km). We 
checked that waves extracted at different radii r > 100 
show the proper fall-off behavior. The waveform is plot- 
ted against the simulated time instead of the more ap- 
propriate retarded time, see below. We present metric 
waveforms only with respect to the retarded time. From 
the plot one clearly identifies the inspiral phase followed 
by the emission related to the HMNS oscillations and 
then the collapse. At early times the well known initial 
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FIG. 21: Metric waveforms from G2P14 evolution. The figure 
shows the real part (top panel) and the imaginary part (bot- 
tom panel) of the r /122 waveform extracted at r = 200 from 
H3 run. The amplitude is also shown in both panels. The 
waveforms are plotted versus the retarded time i re t = t — r* 
where r» is the tortoise radius corresponding to r. The met- 
ric waveform is computed with the FFI and cutoff frequency 
u = 0.002. 



"junk" radiation can be seen. 

Metric waveforms. The "04 waveform is the second 
derivative of the metric waveform h = h+ — 1 h x , which 
represents the actual GW degrees of freedom. The inte- 
gration of the relation 

htm = ipj m (59) 

to obtain h multipoles fr om those of "04, requires some 
attention. In [([l], l67l [o8l Il58l Il59j h is computed via a 
direct (time domain) integration on the simulated time 
domain. The result is affected by a polynomial drift that 
must be corrected by fitting. We refer to this procedure 
as the corrected time domain integration (CTI) . The drift 
observed in the raw integration is not only the linear con- 
tribution expected from the integration of Eq. (|59[) but 
a generic polynomial. It originates from the integration 
of high-frequency noise in the data and has a stochastic 
nature (69|. Recently [H?} an improved CTI procedure, 
which amounts to subtracting a post-Newtonian behav- 
ior before fitting the polynomial correction, which we did 
not consider here, has been developed. In [69[ the fixed- 
frequency integration (FFI) method has been proposed. 
It is based on a spectral integration in the Fourier basis 
and employs a simple high-pass frequency filter against 
spectral leakage. 

We computed /122 with both the CTI and the FFI pro- 
cedure finding comparable results. Differences are dis- 
cussed below. We focus for the moment on I122 com- 
puted with FFI: the integration method does not influ- 
ence the following statements. Fig. [21] shows the met- 
ric waveform r computed from r 1^22 as a function 



of the retarded ti me t re t = t — r * where r* ~ 222 is 
the tortoise radius |l70| corresponding to r = 200. Real 
part, imaginary part and amplitude arc shown for run 
H3. The peak of the amplitude of /122 formally defines 
the merger time, t m , e.g. [13, Il60j . Considering dif- 
ferent resolutions we found t m = 1670, 1710, 1740 and 
1765 (8.23, 8.42, 8.57 and 8.69 ms), respectively, for runs 
H0-3. The corresponding retarded times are im ire t = 
1446, 14861516, 1541 (7.12, 7.32, 7.47, 7.59 ms). The 
metric waveform is composed at early times of six GW 
cycles emitted during the three orbit inspiral. After the 
merger the emission is dominated by the bar-deformed 
HMNS, and the signal has a typical frequency around 
3 kH z wh i ch in creases as the HMNS becomes more com- 
pact [121GI3- Finally, after i ret > 2132 (10.5 ms) the 
GW signal is composed of the quasi-normal-mode ring- 
ing of the BH. We observe from the waves the funda- 
mental frequency ^qnm ^6.5 kHz, compatible with the 
estimate of the BH mass and spin from the apparent 
horizon [13, [H2, i.e. i/ QN m ~ 3.23(10/M BH ) [1-0.63(1 - 
obh) ' 3 ]- The value ^qnm can be estimated by fitting the 
plateau of the GW frequency when it reaches its absolute 
maximum (see Fig. l24[) . A better estimate is provided by 
the frequency of the 1^4 waveforms because the signal is 
less noisy and not contaminated by the integration pro- 
cedure. 

We comment now on the two integration methods em- 
ployed in the post-processing of ^4 to obtain h. In 
the FFI we use as cut off frequency the value vq — 
0.002 (406 Hz), slightly below the GW frequency of 
the initial data. The polynomial correction employed 
in the CTI is a 3rd order polynomial. This choice is 
preferred against a linear or quadratic correction since 
it minimizes experimentally the drift in the raw inte- 
grated waveforms and the oscillations in the modulus 
(see below). The phase difference between the two met- 
ric waveforms amounts to A<f> < 0.06 from early times 
until the collapse. This value is very small and can be 
considered insignificant (see error estimates below and 
Tab. |V|. On the other hand the amplitudes do differ 
in a relevant way. The relative difference is less than 
AA/A < 5% during the inspiral and grows to about 15% 
before the collapse. While these numbers do not seem 
dramatic another important effect is observed. Fig. [22] 
shows the amplitude of r I122 computed with the two 
methods: solid red line for FFI and dashed blue line for 
CTI. The amplitude computed with the CTI shows os- 
cillations during the inspiral. These oscillations converge 
away considering higher resolutions but we were not able 
to remove them completely (the oscillation amplitude is 
larger in waveforms from lower resolutions runs). Waves 
extracted at larger extraction radii show smaller oscilla- 
tions. The use of a linear polynomial correction results in 
very large oscillations (black dotted line) and a prominent 
drift in the waveforms. By contrast the amplitude ob- 
tained by the FFI waveform integrated with i^o = 0.002 is 
free from these oscillations for all the resolutions consid- 
ered. When the FFI with a much lower cutoff frequency, 
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FIG. 22: Gravitational wave amplitude |r/i22| from G2P14 
evolution. The figure shows the amplitude computed with 
the FFI u = 0.002 (red solid line) and u = 0.00035 (green 
dashed-dotted line) and with the CTI where a cubic (blue 
dashed line) and linear polynomial (black dotted line) correc- 
tion is used. 



vq = 0.00035 (71 Hz), is employed the oscillations appear 
again (green dashed-dotted line). Note that a frequency 
of 71 Hz roughly corresponds to the finite length of the 
signal. As observed in [67l.l69l] the oscillations are an un- 
physical effect not related to eccentricity. We mention 
however that we find here a correlation between the os- 
cillations in the amplitude of r /122 and those seen in the 
coordinate separation (Fia ll8[) in the runs with different 
values 77. A gauge effect on the extraction spheres ampli- 
fied in the calculation of h is thus a possible explanation. 
A proper choice of v§ in the FFI can mostly eliminate 
them, while polynomial fitting in CTI is not as robust 
and not performing as well for our data. 

Thermal effects. We discuss now differences between 
the waveforms produced in G2P14 and G2P14hot evo- 
lution. The results refer to HI runs and they are dis- 
cussed referring to Fig. [531 As discussed in Sec. IVIBI 
the simulations show differences of numerical origin al- 
ready in the inspiral. The phase difference between the 
/122 waveforms increases monotonically during the evo- 
lution until t TCt ~ 12 ms when the emission of G2P14 
is practically zero. At the time the two stars experi- 
ence the first contact the G2P14hot /122 waveform has 
accumulated a dephasing of A</> = +0.67 rad and the 
amplitude is about 6 % smaller. Since these differences 
are compatible with the truncation errors for this reso- 
lution and they are not expected in the continuum limit 
(the inspiral is an isentropic process from the fluid point 
of view) they likely have a numerical origin. After the 
contact, thermal effects drive a very different dynamics 
compared to the isentropic case as discussed in Sec. lVIBl 
(see Fig. ITS")) . The dephasing reaches A(p = +2.11 rad at 
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FIG. 23: Comparison between waveforms from G2P14 and 
G2P14hot evolutions. The top panel shows the real part of 
/122 extracted at r = 200 for the two evolutions. The bottom 
panel shows the amplitudes. Vertical lines mark the time 
of the first contact between the stars and the merger time 
for G2P14 evolution. The waveforms are plotted against the 
retarded time t rct =t — r*. Data refer to run HI. 



the merger time and the amplitude is about 14 % smaller. 
Similarly for the r -0| 2 waveform we found a dephasing of 
Acj) = +0.67 rad and a factor —15 % in amplitude at the 
contact time and a dephasing of Acf> = +2.43 rad and a 
factor —50 % in amplitude at the merger time. The post- 
merger emission in the G2P14hot evolution is dominated 
by the emission of the bar-deformed HMNS. A lower fre- 
quency modulation of the signal is visible in the ampli- 
tude of the waveform and corresponds to the nonlinear 
quasi-radial pulsations shown in Fig. [T3] and previously 
discussed. The high-frequency part of the GW signal is 
instead dominated by m = 2 non-axisymmetric nonlin- 
ear modes and it is composed of few frequencies around 
2.7 kHz. A Four ier a nalysis shows results qualitatively 
compatible with |l63j . 

A relevant quantity directly connected to the dy- 
namics of the system is the GW frequency defined as 
W22 = — S(ft, 2 2/ft.22), i-e. the derivative of the GW phase. 
Fig. [53] displays the dimensionless M W22 as computed 
from G2P14 and G2P14hot evolutions. During the in- 
spiral the GW frequency increases monotonically from 
M UJ22 ~ 0.056 to Mw 22 ~ 0.126. This value (green hor- 
izontal solid line, see the bottom panel) is common to 
both evolutions but occurs at different times due to the 
accumulated phase. After the merger, a first local maxi- 
mum is observable with again very close values for G2P14 
(M0J22 ~ 0.181) and G2P14hot (M w 22 ~ 0.184) and 
different times (£ re t ~ 7.84 ms for G2P14hot and t YOt ~ 
8.68 ms for G2P14hot). The GW frequencies present the 
same behavior until this point, thermal effects generate 
only a time shift (retardation). After the merger, the GW 
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frequency of G2P14 presents one oscillation and increases 
from M«22 ~ 0.2 to M 1022 ~ 0.245 when the apparent 
horizon forms. With a steep gradient it further increases 
to MW22 ~ 0.57 which correspond to the fundamental 
QNM of the Kerr BH (i/qnm ~ 6.5 kHz). In the case of 
G2P14hot the GW frequency reflects the dynamics of the 
HMNS: it increases almost linearly with large oscillations 
corresponding to the HMNS quasi-radial oscillations. Af- 
ter the local minimum of the last oscillation the collapse 
happens at a frequency M LO22 ~ 0.27. The QNM fre- 
quency on the G2P14hot evolution is slightly below the 
corresponding iscntropic evolution, ^qnm ~ 6.45 kHz. 
Interestingly the GW frequency becomes negative after 
the first local maximum after the merger in both evo- 
lutions and one oscillation later in G2P14hot (after the 
formation of the apparent horizon in the G2P14 evolu- 
tion). It is difficult to say if these features have a physical 
or a numerical origin. However, the GW frequency zeros 
and negative spikes, roughly correspond to times at which 
the distribution of the rest-mass density in the equato- 
rial plane is almost-spherical. Note in addition that the 
amplitude of the GWs drops to zero (bottom panel of 
Fig. [23]). 

Convergence. We discuss now convergence of the 
waveforms produced in G2P14 evolutions and runs HI, 
H2 and H3. The real part of r^| 2 is reported in Fig. [25] 
for the different runs. The convergence study is focused 
on the inspiral part of the wave. Similar results are found 
for the series HO, HI and H2, with larger absolute errors. 

Fig. [26] displays the (logarithmic) differences between 
the r 1P22 amplitude (top panel) and phase (bottom 
panel) . The difference between H2 and H3 scaled for sec- 
ond order convergence is also plotted. The vertical line in 
the figure marks the merger time as computed from the 
waveforms extrapolated in resolution. The initial junk 
radiation is also cut out from the figure. We observe a 
quite clear 2nd order convergence in the amplitude while 
for the phase the convergence appears slower. The phase 
error is also the dominating error in the waveforms at 
different resolutions. A direct inspection of Fig. [25] al- 
ready highlights this fact. The experimental convergence 
rate measured for the phase is a factor between 1 and 2, 
which suggests that the simulations have "just entered" 
a convergent regime but more resolution or higher order 
numerical methods are required. A proper time shift of 
the waveforms may locally eliminate the phase error thus 
improving the results. We prefer not to perform such 
procedure in order to keep the analysis simple and clean. 
The practice of aligning waveforms for convergence tests 
has been abandoned in some BBH simulations even in 
the more co mpli cated cases of unequal masses and spins, 
see e.g. P. UM flB^ . 

Convergence in the waveforms is lost soon after the 
merger. In the post-merger phase a first order in 
norm convergence is observed in several global quanti- 
ties, cf. Fig. [T?] and the discussion in the previous sec- 
tion. From our data it is not possible to estimate a pre- 
cise point- wise convergence rate in the waveforms: the 
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FIG. 24: Gravitational wave frequency from G2P14 and 
G2P14hot evolutions. The figure shows M U122 computed from 
/122 waveforms. In the top panel the blue solid line refers to 
G2P14 and the red dashed line to G2P14hot. The vertical 
black solid lines mark respectively the merger time and the 
apparent horizon formation for G2P14. The vertical black 
dashed lines mark respectively the merger time and the ap- 
parent horizon formation for G2P14hot. The bottom panel 
gives a detail of the top panel showing the inspiral part. The 
horizontal green solid line marks the value of M U22 at the 
merger for G2P14. Data refer to run HI. 



convergence in amplitude smoothly drops down to first 
order, the phase shows over-convergence. We observe 
however a monotonic dependence on resolution as ev- 
ident from Fig. [25] The reason for this behaviour is 
again the truncation error of the HRSC scheme in this 
strongly nonlinear hydrodynamical phase (characterized 
by a rapid variability in space and time) which dramati- 
cally affects the waveforms. 

Using the Richardson method we extrapolate the in- 
spiral waveforms in order to estimate errors in amplitude 
and phase. Results concerning the maximum error esti- 
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FIG. 25: tf)4 waveforms from G2P14 evolution at different res- 
olutions. The figure shows the real part of the r ip22 waveform 
extracted at r = 200 from runs Hl-3. 



mated during the inspiral are reported in Tab. |V] Using 
the four resolutions and assuming 2nd order convergence, 
we obtain ip4 waveforms with a maximum phase error 
of max (50 ~ 0.3 rad and m&xSA/A ~ 7 %. Assuming 
1st order convergence gives instead max 5(f) ~ 1 rad and 
max 5 Aj A ~ 24 %. Similarly using only the three highest 
resolutions gives max 5(f) ~ 0.6 rad and max 6 A/ A ~ 14 % 
for 2nd order and max<50 ~ 1.2 rad but an unacceptable 
amplitude error for the 1st order assumption. The errors 
on the metric waveform are comparable while generically 
smaller. Furthermore we note that the assumption of 
1st order convergence leads to an evidently non realis- 
tic estimate of the merger time from the extrapolated 
data t m = 2280 (11.23 ms) (Hl-3 data) while it gives a 
more reasonable t m = 1800 (8.87 ms) for 2nd order (Hl-3 
data). 

Some care is needed in interpreting these results be- 
cause of the previous discussion. However, from our re- 
sults we conclude that: (i) the series HO-2 while showing 
convergence is too inaccurate and not reliable for error 
estimates; (ii) assuming 1st order convergence for the se- 
ries Hl-3 and HO-3 is not appropriate and overestimates 
the actual errors; (iii) a conservative and realistic error 
estimate is provided by Hl-3 data assuming 2nd order 
convergence; (iv) the error on the extrapolated data can 
be estimated as the difference between the extrapolation 
with four resolutions (HO-3) and three resolutions (Hl-3), 
e.g. max 5(f) ~ 0.33 rad and m&x5A/A ~ 7.4 % for r^| 2 . 
As indicated in [53], the main error on the waveforms is 
due, as expected, to resolution rather than finite radius 
extraction. Since our results are clearly dominated by 
truncation error we did not investigate finite extraction 
errors but left that study for future work. 
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FIG. 26: Convergence of phase and amplitude in rip22 wave- 
form from G2P14 evolution. The figure shows the difference 
between runs HI and H2 (blue solid line), between runs H2 
and H3 (green solid-dotted line) as well as the difference be- 
tween H2 and H3 (red dashed line) scaled for 2nd order con- 
vergence. Top panel displays the logarithm of amplitude dif- 
ferences, bottom panel displays the logarithm of phase differ- 
ences. The vertical line indicates the merger. 



VII. CONCLUSIONS 

In this work we presented detailed tests and first full 
scale evolutions for a new computer code aimed at 3+1 
numerical studies of general relativistic matter space- 
times. The implementation represents an upgrade of the 
BAM code developed previously for vacuum systems (62T - 
[63 | . The code solves the flux-conservative formulation of 
the GRHD equations [65[ coupled with the BSSNOK sys- 
tem. Mesh refinement (moving boxes) and a metric solver 
were already provided by the BAM code. The GRHD 
equations are solved by means of robust HRSC meth- 
ods [94|, [96H99I ] and share the same grid and the same 
time stepping algorithm with the metric solver. We de- 
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TABLE V: Error estimates during inspiral for extrapolated 
waveforms. Columns: data used for extrapolation, assumed 
order of convergence, waveform, maximum absolute error in 
phase, maximum relative error in phase, maximum relative 
error in amplitude. 



data 


r 


waveform 


max 8(f) [rad] 


max 5(j)/(j) 


'%] max 5 A/ A [%] 




2 


r^22 


0.29 


0.68 


6.89 


HO-3 




rh 2 2 


0.26 


0.64 


2.65 




1 


r^22 


1.04 


2.40 


24.02 




rh 2 2 


0.87 


2.17 


13.31 




2 


r^22 


0.62 


1.43 


14.30 


Hl-3 




rh 22 


0.53 


1.32 


6.85 




1 


r^22 


1.20 


2.85 


> 100 




rh 2 2 


1.10 


2.81 


^ 100 



are probably close but not still not optimal for produc- 
tion runs, we found 2nd order convergence during the 
inspiral phase without any time-shifting procedure. For 
the first time in BNS simulations we consistently esti- 
mated error-bars on the waveforms by extrapolating the 
results in resolution. 

This work represents our first contribution to the study 
of GWs from BNS systems. The methodology presented 
here will be applied in future works. In particular, we 
plan in the near future the production of accurate and 
convergent inspiral waveforms from initial configurations 
described by different realistic EoS in order to investigate 
the impact of the EoS on the GW signal @, [HI, HE H3] • 



scribed in detail the equations and the numerical method 
implemented. The code allows the use of one-parameter 
EoS tables and implements the hybridization procedure 
of 36] to model thermal effects. We proposed a sim- 
ple, thermodynamically consistent interpolation scheme 
for one-parameter EoSs. 

We validated the code against a number of stringent 
tests involving single star spacetimcs. We studied the 
performance of different reconstruction procedures and 
the convergence of the code. We found that in our set up 
the use of CENO reconstruction shows a clear 2nd order 
convergence and leads to smaller global truncation errors 
compared to the other methods. 

We presented test evolutions of irrotational equal-mass 
binary neutron star configurations. Matter is described 
by simple polytropic and ideal gas EoSs. In both evolu- 
tions the formation of an HMNS is observed. In the isen- 
tropic case it collapses quite rapidly producing a Kerr BH 
surrounded by a disk of mass Md < 2 x 10~ 2 Afo while 
in the other case the HMNS survives for 9 ms due to 
thermal pressure support. 

We investigated the impact of different reconstruction 
methods on the dynamics of the HMNS. Our results high- 
light a strong influence of the numerical method and of 
the resolution on the simulated physics. Precise quantita- 
tive statements on the post-merger phase require extreme 
care. 

We presented new results concerning the use of dif- 
ferent values of the damping parameter r/ in the shift 
condition. Small values of r\ are found to produce a less 
coordinate-related eccentricity during inspiral and to re- 
duce the coordinate size of the final BH. 

The gravitational radiation emitted by the system 
was investigated. We characterized the GWs (metric 
waveforms) and discussed two different methods to com- 
pute them from ip4 waveform: a time domain (CTI, 
e -g- USUI]) an d Fourier domain (FFI, integration. 
The FFI provides better results minimizing unphysicaJ 
drifts and oscillations. We performed self-convergence 
tests of the waveforms. While the resolutions employed 



Acknowledgments 

The authors sincerely thank David Hilditch and 
Alcssandro Nagar for many valuable discussions. The 
authors thank Dorccn Mullcr for discussion on metric 
waveform integration. The authors thank the Mcudon 
group for making publicly available L0RENE initial data 
and Eric Gourgoulhon for explanations. During the de- 
velopment of the BAM matter code and of this paper the 
authors benefitted from valuable discussions and corre- 
spondence with Luca Baiotti, Thibault Damour, Roman 
Gold, Roberto Dc Pietri, Harald Dimmclmcicr, Luca Del 
Zanna, Frank Lofficr, Pedro Montero, Luciano Rezzolla, 
and Wolfgang Tichy. This work was supported in part by 
DFG grant SFB/Transregio 7 "Gravitational Wave As- 
tronomy" . Computations where performed at LRZ (Mu- 
nich) and on local clusters at TPI (Jena). 



Appendix A: Algorithm to recover primitives 

In this appendix we summarize the procedure and 
the equations to recover the primitive variables from 
the conservative ones. The specific algorithm adopted 
has been deve l oped in a number of previous publica- 
tions HI [Tol, Gjjl fl20l fl36l fl65ll . Other algorithms 
have been developed, e.g. |l66l . ll67j? in particular we have 
successfully tested the method of jl66| in single star evo- 
lutions. We leave for future work extensive tests and 
comparisons between different algorithms. 

We first discuss an iterative proc edur e for a general 
EoS, p = P{p, e), first employed in 165]. Specific pro- 
cedures can be designed once a specific form of EoS is 
given [92]]. We then present the procedure we adopt in 
the case of c old E oSs which is based on an iterative algo- 
rithm for p [llOj . Finally we describe the modifications 
introduced to handle the presence of the artificial atmo- 
sphere. 
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Following [1061 ] we invert Eqs. (f2"2"j) to find 
5' 



the latter is determined by 



v\p) 
W{p) 



T + D+p ' 

T+p + D 



pip) 

e(p) = Z^ 1 



(r+p + D) -S 2 



(r+p + Z?) 2 -^ 2 - Wp-D 



(Al) 
(A2) 

(A3) 
(A4) 



Primitive variables can be computed from the conserva- 
tive variables using the above relations once the pressure 
is known. The pressure is determined by the EoS looking 
for the root of the nonlinear algebraic equation 

f{p)=p-P(p{p),e(p)) = . (A5) 

The algorithm used is the Ncwton-Raphson iteration, 

old fip) 



p ncw = p° 1 ' 



fiP) 



The derivative of / is given by 

f{p) = l - X d-p- K dp-> 
DS 2 



dp 

dp ~ (D + p + t) 2 ^{D + p + r) 2 

de pS 2 

dp 



D{(D +p + r) 2 - S 2 ) 



3/2 



(A6) 

(A7) 
(A8) 
(A9) 



In the case of a one-parameter EoS, p = P(p), h = h(p), 
e = e(p), and 



W 



S 2 



(Dhf 



(A10) 



are functions of the density and the conservative variables 
only. Primitive variables can be computed from p, once 



g(p) = W(p)p - D , 



using again a Newton-Raphson root finder with 



(All) 



g'(p) 
h'(p) 



W(p) - p 



S 2 h'{ P ) 
WD 2 h? 



P P 



(A12) 
(A13) 



Note that h'(p) = ^ if the principle of thermodynamics 
at zero temperature is applied. 

The recovery procedure for a general EoS can fail at 
low densities, e.g. in presence of the atmosphere. A rea- 
son for this is simply machine accuracy: since typically 
p -^i D, for very low densities the Newton-Raphson algo- 
rithm docs not converge. A further complication is the 
spuriously high value of the velocity generated by the ar- 
tifical atmospher e tre atment. To handle these problems, 
similarly to |H, [TToj . we combine both the algorithms 
described above and we implement a set of hierarchic 
prescriptions which are found to work in practice. Specif- 
ically: (i) every time a grid point reaches a density below 
the atmosphere threshold density the code sets atmo- 
sphere values both in primitives and conservatives and 
continues to the next point; (ii) if the general EoS algo- 
rithm does not converge due to a too small value of the 
pressure, i.e. p ncw — p° ld < e ~ 10~ 12 , then atmosphere 
values are set; (iii) if it returns unphysical values for e, p 
or p then the code tries the inversion with the algorithm 
for the cold EoS; (iv) if it returns unphysical values of 
v 2 then atmosphere values are set; (v) if the algorithm 
for the cold EoS does not converge or returns unphysical 
values not cured by finer grid levels, then the code stops. 
In practice we observe the necessity of (ii), (iii) and (iv) 
only in binary simulations. 
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